Assessment of the Impact of GNSS Processing Strategies on the Long-Term Parameters of 20 Years IWV Time Series

: Advanced processing of collected global navigation satellite systems (GNSS) observations allows for the estimation of zenith tropospheric delay (ZTD), which in turn can be converted to the integrated water vapour (IWV). The proper estimation of GNSS IWV can be affected by the adopted GNSS processing strategy. To verify which of its elements cause deterioration and which improve the estimated GNSS IWV, we conducted eight reprocessings of 20 years of GPS observations (01.1996–12.2015). In each of them, we applied a different mapping function, the zenith hydrostatic delay (ZHD) a priori value, the cut-off angle, software, and the positioning method. Obtained in such a way, the ZTD time series were converted to the IWV using the meteorological parameters sourced from the ERA-Interim. Then, based on them, the long-term parameters were estimated and compared to those obtained from the IWV derived from the radio sounding (RS) observations. In this paper, we analyzed long-term parameters such as IWV mean values, linear trends, and amplitudes of annual and semiannual oscillations. A comparative analysis showed, inter alia , that in terms of the investigation of the IWV linear trend the precise point positioning (PPP) method is characterized by higher accuracy than the differential one. It was also found that using the GPT2 model and the higher elevation mask brings beneﬁts to the GNSS IWV linear trend estimation.


Introduction
Global navigation satellite systems (GNSS) find many uses in the field of navigation, geodesy, and geophysical research. With regard to geophysical studies, more and more effort is being put into developing methods related to climate monitoring. The usage of GNSS in atmosphere studies results from the fact that the advanced processing of GPS observations allows for estimating zenith tropospheric delay (ZTD), which can be potentially used for climate monitoring [1,2], in water vapour content estimation [3,4] or in meteorology applications [5]. This parameter contains the impact of both the hydrostatic (zenith hydrostatic delay, ZHD) and the wet (zenith wet delay, ZWD) parts of the troposphere on the GPS radio signal propagation. ZHD is characterized by slight changes in time and space, and therefore is relatively easy to model. The opposite situation occurs in the case of ZWD, because its time and space variation exhibits a similarity with the water vapour content, which is one of the most important natural greenhouse gases in the atmosphere [6]. Consequently, ZWD can be converted into integrated water vapour (IWV) or precipitable water vapour (PWV) [7,8]. Thanks to the high temporal and spatial resolution of GNSS IWV, this parameter is used in many studies related to the climate [9,10], weather forecast [11], or very long baseline interferometry of VLBI [12]. According to Ning et al. (2016) [13], the uncertainties in ZTD dominate the error budget of the IWV at the level of 75%. Therefore, it is extremely important to determine the ZTD value as reliable, because its accuracy directly affects the quality of IWV. From the beginning of GNSS permanent measurements, various methods and models have been applied to observations processing. Each of them can affect final ZTD solutions in a different way and, in consequence, the GNSS IWV/PWV time series.
ZTD is one of many unknowns that are estimated during the GNSS positioning. Its value is estimated in several steps and the first of them is adopting ZHD a priori value. Subsequently, the adopted value is mapped to the satellite direction through a mapping function dedicated to the hydrostatic delay. Then, the remaining part of the delay towards the satellite is mapped to the zenith direction using the mapping function dedicated to the wet delay. In the final step, the correction of the firstly adopted ZHD value is estimated. If the firstly adopted ZHD a priori value is correct, then the obtained correction can be identified with ZWD. The sum of ZHD and ZWD gives ZTD. Since they are various mapping functions and sources of ZHD a priori value, many studies were focused on the analysis of their possible influence on ZTD estimation. Tregoning and Herring (2006) [14] have compared ZHD from the global pressure and temperature model (GPT) [15], the standard temperature and pressure model (STP), and in the situ measurements near the GPS station. At the same time, Vey et al. (2006) [16] tested using the empirical Niell mapping function (NMF) [17] and, based on the numerical weather model (NWM), isobaric mapping function (IMF) [18]. Their work was undertaken inter alia by Steingenberger et al. (2009) [19] who have analyzed the impact of the commonly used global mapping function (GMF) [20] and Vienna Mapping Function 1 (VMF1) [21], with sourcing data from GPT and the European Centre for Medium-Range Weather Forecasts (ECMWF). In all those cases, the differences in the ZTD time series obtained according to various processing schemes were found. The fact that these kinds of differences may translate into GPS IWV time series was confirmed by Thomas et al. (2011) [22], who, next to the ZHD and the mapping function (NMF/STP compared to VMF1/ECMWF), have also verified the influence of applying the atmospheric loading model (ATML) and individual antenna calibration. The biases between the ZTD time series obtained according to the above strategies caused differences in GPS PWV. The influence of elevation-dependent errors such as antenna phase center variations [23] or the signal multipath [24] on the GPS PWV was in turn undertaken by Vey et al. (2009) [25], who have verified that changing the elevation mask from 3 • to 15 • caused 0.75 mm differences in obtained PWV values.
All these studies were based on up to several years of observations. Therefore, it was hard to assess whether by investigating them differences may also arise in case of long-term changes, or not. The fact that various processing strategies can result in various values of the ZTD linear trend was presented by Baldysz et al. (2016) [2]. Based on the reprocessing of 18 years of observations according to two different processing strategies, they have shown that the differences in the obtained linear trends reached up to 0.26 mm/year. Due to the lot of differences between these two strategies, it was impossible to assess which factor was most responsible for the observed changes. Studies on this issue have been taken by Dousa et al. (2017), who compared seven different reprocessing variants [26]. They differ from each other in terms of mapping function, cut-off angle, non-tidal atmospheric loading, higher-order ionospheric corrections, and sampling troposphere gradients. On the basis of 172 EPN (EUREF Permanent Network) [27] stations and at least ten years of data (max span was 1996-2014), they have verified, e.g., that ZTD linear trend values may vary from each other up to 0.12 mm/year in the case of a different mapping function and up to 0.35 mm/year in the case of using various elevation masks. These results were not compared with other in situ measurements. Therefore, it is not known which of them are more reliable. This was partly verified by Ning et al. (2016) [28], who, on the basis of 20 years of GPS observations, have shown that applying various elevation masks during GPS observations processing causes differences in long-time changes of GPS IWV. They verified that higher correlations between GPS IWV and IWV from ERA-Interim [29] occurs when a 25 • elevation mask is used instead of a 10 • . However, it should be kept in mind that ERA-Interim is a numerical model, and although it contains reprocessed and global coverage data, there is still lack of comparison with the in situ meteorological measurements.
In order to verify the influence of other elements of the GNSS processing strategy on the long-time changes of the GNSS IWV time series, we conducted eight reprocessings for the selected EPN stations. Among them were a mapping function, an ZHD a priori value, an elevation mask, software, and the positioning method, all of which were taken into account. Using Bernese GNSS Software 5.2 [30] and GAMIT 10.50 [31], we tested 20-years ZTD time series (01.1996-12.2015) from the 20 EPN stations. As a troposphere model, VMF1, GMF/GPT, GPT2, and NMF/STP, described in detail in section two, were used. The ZHD and GNSS IWV values were obtained on the basis of pressure and mean temperature (T m ) data from the ERA-Interim. The validation of this model with in situ measurements, as well as the estimation of errors which it can translate to GNSS IWV, are also given in section two. The final results, including mean values, annual and semi-annual oscillation, and linear trends, are presented in section three, which also contains a comparison of the observations from the radio sounding (RS) stations.

Atmospheric Propagation Delay
Electromagnetic wave delay (T e ) is caused by refraction and attenuation in the troposphere and is defined by the following formula [7]: in which c is the speed of light in the vacuum, τ is the delay measured in the unit of time, and N is the neutral atmospheric refractivity, which can be divided in hydrostatic (dry gases) and wet (water vapour) components [32]: in which P is the total air pressure in [hPa], T is the temperature in [K], e is the water vapour pressure in [hPa], and K 1 , K 2 , K 3 are the air refractivity parameters. In our study, we used their values estimated by Rüger (2002) [33], which amount to 77.689 ± 0.0094 [1/hPa], 71.295 ± 1.3 [K/hPa], and 3.75463 × 10 −5 ± 0.0076 × 10 −5 [K 2 /hPa] respectively. In many related studies, the ZTD parameter (often used interchangeably with "zenith path delay", abbreviated as ZPD) is analyzed as: Both parameters can be given by the following integral along the atmospheric profile: in which in which ϕ is the ellipsoidal latitude, and h is the height in [m]. IWV can be obtained using integral along the atmospheric profile: or using ZWD: in which R v is the specific gas constant of water vapour and amounts to 461.5 [J/kg K]. Π(T m ) is a dimensionless quantity that can be obtained from the mean weighted temperature of the atmosphere T m : Several empirical relations are given to determine T m using linear relationships with surface temperature (T s ), which are tuned to the specific area and season [7,[35][36][37][38]. This method is easy to use but does not provide high accuracy, e.g., Bevis et. al. (1994) [35] obtained the T m − T s RMS error at the level of 4.74 K. The relative error of the parameter Π closely approximates the relative error of T m . The 1 K error of T m caused 0.41-0.32% error of Π for T m value between 240-300 K. This value translates into the IWV error, which amounts to 0.05-0.1 kg m −2 for the ZWD value between 0.1-0.2 m. To obtain high accuracy of T m , the following formula can be applied: This formula can be only used when profiles through the atmosphere in zenith direction are available. So, it is necessary to have the RS observations or data from the numerical weather model (NWM), which in many cases may not be available, especially in the case of the RS stations that have low spatial and temporal resolution. IWV contains the same information as PWV. The only difference between these two parameters results from the fact that PWV is adjusted by water density and therefore is expressed in the unit of mm.
Both Equations (8) and (9) give the same results (if T m is obtained using Equation (11)), but the second one can be applied using ZWD estimated from GNSS processing. Thanks to this, GNSS data can be used to study the variability of the water vapour content in the atmosphere.

ZTD from GNSS Processing
The total tropospheric delay in direction of the satellite-slant tropospheric delay (STD) can be divided into hydrostatic (SHD) and wet (SWD) components. Both of them can be determined using mapping functions and the corresponding delays at the zenith using the following formula: in which SHD and SWD are slant total delays (in the direction to the satellite) for the hydrostatic and wet part of the tropospheric delay, which are projected into zenith direction and vice versa using hydrostatic and wet mapping functions, m f h and m f w , which depend on the satellite elevation angle (el) and can be described as follows [39]: 1+c sin(el) + a sin(el)+ b sin(el)+c (13) Remote Sens. 2018, 10, 496 5 of 25 in which a h , b h , c h , a w , b w , c w are coefficients for hydrostatic and wet delays. Nowadays, mapping functions based on the concept of the Vienna Mapping Function or Global Mapping Function are commonly used in GNSS data processing. G N and G E are azimuth (α)-dependent gradients expressed in north and east direction, and m f g is the gradient mapping function. The Chen and Herring gradient mapping function [40] is mostly used in GNSS data processing.
In our study, we used GPS observations from 20 EPN stations ( Figure 1). Their choice was motivated by (i) an as long as possible time span of collected observations and (ii) the availability of an RS station in the range of 150 km. Selected stations operated constantly from 1996 until the end of 2015, therefore delivering 20-years of observations, which we processed according to eight different processing strategies. Four of them were conducted in Bernese GNSS Software 5.2 (with 'BSW' acronym) and were focused on the troposphere modelling and positioning method. For this purpose, two commonly used mapping functions with troposphere models were applied. The first of them was routinely used in EPN combination of VMF1 and ECMWF and the second was GMF/GPT. These two approaches differ from each other in terms of both mapping function and source of ZHD a priori value. In the case of VMF1 b w and c w , coefficients were adopted from NMF, and b h , c h were calculated on the basis of one year of data from ERA-40 [41] through the ray-tracing procedure at 3.3 • initial elevation angle (b h assumed the constant value and c h was described by the model related to day of year and latitude). The last two coefficients of VMF1 a h , a w were delivered through a ray-tracing procedure of operational weather forecasts of ECMWF at 0, 6, 12, and 18 UTC. In GMF, b h , c h , b w and c w were adopted from VMF1, whereas a h and a w were determined by a single ray-tracing of ERA-40 at 3.3 • initial elevation angle for the 1999-2002 period of time. Therefore, the only difference between these two mapping functions occurred for a h and a w coefficients, which in VMF1 are routinely calculated and in the GMF are modelled. The same difference occurred for the ZHD a priori value. In the first combination (VMF1/ECMWF), it is sourced from the operational forecast of ECMWF, whereas in the second combination (GMF/GPT) it is estimated from the empirical model GPT, which is constructed from the same set of data as a h and a w in the GMF. These two approaches in troposphere modelling were applied in both the precise point positioning (PPP) and double difference (DD) method, and the solutions are named BSW_PPP_VMF, BSW_PPP_GMF, BSW_DD_VMF, and BSW_DD_GMF (see details in Table 1). PPP is an absolute method that allows for direct estimation of station position, whereas DD is a relative approach in which baselines between stations are estimated. Using PPP method, each station is analyzed separately, so the accuracy of final solutions are dependent mainly on the local conditions. Thanks to that, the station's position is not affected by network related errors. However, the quality of PPP solutions strongly depends on the quality of satellite ephemerides, as well as on used models. DD method is less sensitive to quality of mentioned products, because a lot of errors are reduced during relative positioning. The main disadvantage of DD is the fact that, in contrast to the PPP, it can be affected by network related errors that can be transferred through baselines. In both these methods, conception of tropospheric delay estimation is the same.
The remaining four strategies were carried out in GAMIT 10.50 software. In the first of them (GAM_DD_VMF) we decided to apply the same combination of VMF1 and ECMWF as was used in BSW_DD_VMF in order to verify whether the internal processing scheme of both software may cause differences. This solution is also compatible with Military University of Technology (MUT) contribution to the EPN Repro2 project [42]. In the second one (GAM_DD_STP), the previously often-used NMF and STP combination was applied. The NMF coefficients for hydrostatic and wet delay were calculated based on US SA66 (US Standard Troposphere 1966) and their values are described as a function of day of the year, latitude, and station height. ZHD a priori value was, in this case, estimated on the basis of standard temperature and pressure and the Saastamoinen formula. This choice was motivated by the fact that this troposphere model was used in the first official reprocessing of EPN and may be responsible for some of the differences founded by Baldysz et al. (2016) [2]. The last reprocessing concerned troposphere modelling based on the GPT2 [43] model (GAM_DD_GPT2). This model does not have to be used together with mapping function, because it delivers both a priori ZHD and a h and a w coefficients. It was built on the same numerical model as GPT, but, in contrast to its predecessor, 10 years of data (2001-2010) were used for this purpose. It is also characterized by higher horizontal and vertical resolution. The fourth reprocessing was exactly the same as GAM_DD_VMF, but instead of the 5 • cut-off angle we adopted 20 • and named this solution GAM_DD_EL20.
All of the above reprocessings were conducted using reprocessed ephemerides from CODE [44]. In each case, also higher order ionospheric corrections were included in the calculation process, since their impact on the GNSS solutions was tested [45,46] and finally recommended for IGS or EPN analysis. As an antenna model, we used the type mean models (IGS08), and calculated solutions were expressed in the IGb08 frame. On the basis of obtained hourly ZTD, we calculated the daily ones using only these values whose formal estimation error was within the 3σ criterion. In next step, they were converted to the IWV by using formulas described in Section 2.1 and Equation (13) from Section 2.2. Finally, we estimated the annual and semi-annual oscillations and linear trends for each GNSS station and solution. This was done according to the following model: in which x(t) is the observed time series, v x t is the linear trend, and A O A A I A etc. are coefficient for the condition of sine and cosines using the method of least squares. The number of concerned sine and cosine conditions was chosen individually for each station by making Lomb-Scargle periodograms [47].

Meteorological Parameters
As we mentioned before, to calculate IWV from ZTD the meteorological parameters such as surface pressure (to estimate ZHD from the Saastamoinen model (7) and therefore ZWD (3)) and T m (to calculate IWV from ZWD (11)) are necessary. Unfortunately, there are not many GNSS stations equipped with meteorological sensors, and, even in the case of the ones that are, the quality of this data and its completeness are questionable. Thus, it is a problem to obtain reliable meteorological data on GNSS station location.
In several studies [3,48,49], data from neighboring synoptic stations are interpolated to the localization of the GPS station and then are used in further analysis. Such an approach allows for obtaining surface pressure and T m , which is calculated from the T m − T s relationship. However, obtained in such way T m is not characterized by very high accuracy and therefore may cause errors during the ZWD to IWV conversion. Data from synoptic stations are not suitable for T m estimation according to Equation (11) due to the fact that they do not provide temperature and humidity profiles (or dew point temperature). Such profiles can be obtained using radiosondes or microwave measurements. Unfortunately, their spatial distribution is not sufficient enough to make it possible to use them in all cases. For example, on the Polish territory they are only 3 sounding stations, whereas GNSS permanent stations are more than 400. Therefore, we decided to use meteorological data sourced from the ERA-Interim (Dee et al., 2011) [29]. This model allowed us obtain data from 1979 to the present with about 80 km and 6 h of spatial and temporal resolution, respectively. Its advantage is also the possibility of getting various profiles up to 0.1 hPa level. What is more, the ERA-Interim data are re-analyzed, so it can be stated that they are continuous.
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 25 obtaining surface pressure and , which is calculated from the − relationship. However, obtained in such way is not characterized by very high accuracy and therefore may cause errors during the ZWD to IWV conversion. Data from synoptic stations are not suitable for estimation according to Equation (11) due to the fact that they do not provide temperature and humidity profiles (or dew point temperature). Such profiles can be obtained using radiosondes or microwave measurements. Unfortunately, their spatial distribution is not sufficient enough to make it possible to use them in all cases. For example, on the Polish territory they are only 3 sounding stations, whereas GNSS permanent stations are more than 400. Therefore, we decided to use meteorological data sourced from the ERA-Interim (Dee et al., 2011) [29]. This model allowed us obtain data from 1979 to the present with about 80 km and 6 h of spatial and temporal resolution, respectively. Its advantage is also the possibility of getting various profiles up to 0.1 hPa level. What is more, the ERA-Interim data are re-analyzed, so it can be stated that they are continuous. Based on the ERA-Interim data, we calculate surface pressure ( ) on the GNSS station locations using the following steps: (i) firstly, four grid nodes nearest to the GNSS station are selected; (ii) at these nodes is converted to sea pressure ( ) using Equation (15); and (iii) then, is interpolated at the GNSS latitude and longitude using bilinear interpolation [50] and converted to using transformed Equation (13).
In the case of based on the ERA-Interim profiles, we estimated using Equation (11) and then performed a bilinear interpolation to obtain its value at the GNSS station location. We did not apply any height corrections.
To validate our methodology and to check ERA-Interim accuracy, we performed a comparison to the in situ measurements based on the soundings stations. We selected 20 sounding stations located as close as possible to the analyzed GNSS stations. In Table 2, we present these stations together with the distance to the nearest GNSS stations. The location of the sounding stations are also shown in Figure 1. Based on the ERA-Interim data, we calculate surface pressure (P) on the GNSS station locations using the following steps: (i) firstly, four grid nodes nearest to the GNSS station are selected; (ii) P at these nodes is converted to sea pressure (P sea ) using Equation (15); and (iii) then, P sea is interpolated at the GNSS latitude and longitude using bilinear interpolation [50] and converted to P using transformed Equation (13).
In the case of T m based on the ERA-Interim profiles, we estimated T m using Equation (11) and then performed a bilinear interpolation to obtain its value at the GNSS station location. We did not apply any height corrections.
To validate our methodology and to check ERA-Interim accuracy, we performed a comparison to the in situ measurements based on the soundings stations. We selected 20 sounding stations located as close as possible to the analyzed GNSS stations. In Table 2, we present these stations together with the distance to the nearest GNSS stations. The location of the sounding stations are also shown in Figure 1. Our comparison was made based on the same 20-year time range for which we analyzed the GNSS data. In Table 3, detailed differences between the measured and modeled parameters are shown. It can be noticed that the maximum pressure bias was obtained for station 04018 and amounted to −1.26 hPa. Based on Equation (7), it can be seen that this value corresponds to about −2.9 mm of ZHD. If we look at the absolute mean value of ∆P, we can notice that the impact of wrong modeled pressure on ZHD is smaller and amounts to about 1.1 mm. Better results were obtained for annual and semiannual amplitudes and linear trend in which the absolute mean differences amounts to 0.10 hPa, 0.05 hPa, and 0.02 hPa/year respectively, which is around 0.2 mm, 0.1 mm, and 0.05 mm/year of ZHD. It should be remembered that these values, which can be called errors caused by the meteorological model, translate directly to ZWD and further to the IWV value.
In the case of ∆T m , we see that the biggest bias, 1.71 K, is for station 08160. For T m = 270 K and ZWD = 0.1 m, this value caused about 0.1 kg m −2 differences of IWV. As we can see, the absolute mean value is much smaller (0.55 K) and hence the impact on IWV is also smaller (about 0.03 kg m −2 ). Based on these considerations, we state that the impact of wrong modeled and interpolated v T m can be neglected, because the obtained absolute mean of v T m (0.01 K) translates to IWV at the level of about 0.0006 kg m −2 /year.
To assess how the modeled and interpolated P and T m affect the IWV, we made a calculation based on Equations (7), (9), and (10). Values from Table 3 (∆P and ∆T m ) were used to estimate errors of ZHD (which translates to ZWD (3)) and Π(T m ). Then, using the mean value of T m and ZWD derived from the sounding observations (T m S and ZWD S , respectively), we calculated the differences of mean IWV (∆IWV), the amplitude of annual (∆A A IWV ) and semiannual (∆A S IWV ) oscillations, and the linear trend (∆v IWV ). All results we show in Table 4.
Based on the results presented in Table 4, we state that the used source of meteorological data (the ERA-Interim), as well as the adopted interpolation methods, give results that can be used for further analysis. In our research we mostly pay attention to long-term parameters such as oscillation amplitudes and trend. As it is seen in Table 4, the used meteorological parameters do not affect them strongly (0.05 kg m −2 , 0.02 kg m −2 , and 0.006 kg m −2 /year, for annual and semiannual oscillations, and linear trend, respectively). As it will be seen later in our paper, the obtained absolute mean differences are below the GNSS IWV statistic parameters estimation error. Slightly worse results were obtained for ∆IWV, in which the biggest value is for station 04018. This is not surprising, because for this station also the biggest ∆P was noticed (Table 3). Even so, this value is still acceptable, the more that absolute mean is only 0.16 kg m −2 . In spite of all, due to large distance between GNSS and the soundings stations, and the height differences between them, in our paper we did not analyze in detail the total value of IWV. We believe that these distances do not significantly affect the differences in trend and average amplitudes. Table 3. Long-term statistics (amplitude of annual (A A ) and semiannual (A S ) oscillation, trend (v), and bias) of surface pressure (P) and weighted mean water vapour temperature (T m ) differences between RS measurements and data derived from ERA-Interim model. Estimation errors: σ∆P ∼ = ±0.03 hPa;  Table 4. Differences of mean IWV (∆IWV), amplitude of annual (∆A A IWV ) and semiannual (∆A S IWV ) oscillations, and linear trend (∆v IWV ) caused by differences of meteorological parameters presented in Table 3.

Results
In this section, we present our analysis of the calculated GNSS IWV time series. The mean differences between GNSS IWV and IWV from RS are present in this section. Subsequently, the interannual variations of GNSS IWV, in terms of annual and semi-annual components, and their comparison to the RS IWV are presented. Finally, we would like to draw the greatest attention to the analysis of long-term differences between all conducted reprocessings, in particular to their validation using RS observations. The comparison between GNSS IWV and RS IWV was made only for days in which data were available of both solutions.

Mean IWV
The selected GNSS stations were characterized by various prevailing humidity conditions. The mean IWV values from all analyzed GNSS IWV solutions were in the range of 12.07 kg m −2 to 22.05 kg m −2 ( Table A1 in Appendix A). These results were consistent with those obtained from RS IWV, which were in the range of 11.46 kg m −2 to 19.51 kg m −2 ( Table A2 in Appendix A). Although both GNSS IWV and RS IWV were characterized by a similar mean value, it can be noticed that the RS results were slightly lower than the GNSS ones. This is visible in Figure 2, in which differences between them are given. Only for two stations (MATE, MOPI), the mean GNSS IWV was clearly lower than the mean IWV from the corresponding RS stations. For the MATE and 16320, the height difference is at the level of 480 m, whereas for the MOPI and 11035 it is 335 m. Therefore, the characteristic decrease in air humidity with increasing height can be noticed here. Generally, absolute mean differences between GNSS IWV and RS IWV were smaller than 1.50 kg m −2 (only in the case of one solution it was 1.51 kg m −2 ) ( Table A3 in Appendix A). Among all stations, the biggest similarity between the GNSS and RS results occurred for the GRAZ and 11,240 stations. In this case, all of the Bernese solutions were below 0.25 kg m −2 , and only two GAMIT solutions were higher than 0.50 kg m −2 . This consistency results mostly from the small distance between the GNSS and RS station, which is only 9 km. Analogical situation occurred in the case of HERS and 03882, which had differences in the mean IWV value at the level of 0.57 kg m −2 . These two stations are only 4 km away from each other. On the other hand, PENC, SFER, and BOR1 with corresponding values of 12843, 08495, and 12425 are several dozen kilometers from each other (up to 130 km) and are still characterized by one of the smallest discrepancies with respect to (w.r.t.) RS. The biggest differences occurred for the EBRE and 08160 stations and reached up to 5.37 kg m −2 . These stations are in the 147 km distance at the various station heights (∆h = 200 m). This probably caused reducing of the Mediterranean water masses influence on the local humidity conditions for the inland localized 08160 station. There were also other pairs of GNSS and RS stations that were a large distance from each other, like in the case of BOR1 with 12425 and WARE with 10410. However, in contrast to the EBRE and 08160 stations, they are characterized by much smaller height differences, which were at the level of 27 m and 11 m (Table 3). An interesting situation occurred for MEDI and corresponding 16144. These two stations are 15 km from each other, and the ∆H for them is about 1 m, but still, differences in the mean IWV between them were higher than in the case of the stations separated by a much larger distance. This situation was confirmed by each of the GNSS solutions, so probably there occurred some local conditions that affected the GNSS observations in that way. On the other hand, it has to be also kept in mind that all differences between RS and GNSS may be partly caused by the model uncertainties. As it was shown in Table 4, differences of the mean IWV value between the ERA-Interim and RS observations varied from −0.44 kg m −2 to 0.43 kg m −2 (observations-ERA-Interim). Due to the fact that in our conversion from ZTD to IWV we sourced meteorological parameters from the ERA-Interim, the uncertainties of the model may translate to the uncertainties of those calculated from GNSS IWV. differences of the mean IWV value between the ERA-Interim and RS observations varied from −0.44 kg m −2 to 0.43 kg m −2 (observations-ERA-Interim). Due to the fact that in our conversion from ZTD to IWV we sourced meteorological parameters from the ERA-Interim, the uncertainties of the model may translate to the uncertainties of those calculated from GNSS IWV.

IWV Amplitudes
Because they took place in the mean IWV, the values of the annual oscillations did not significantly vary between the particular GNSS solutions (Figure 3) (Table A4 in Appendix). More precisely, the biggest discrepancies were found for the MOPI station and amounted to 0.33 kg m −2 (GAM_DD_STP and GAM_DD_EL20) ( Table A5 in Appendix). This value makes up about 4% of the total value of the annual component for this station and is the greatest deviation from all solutions both in terms of absolute and relative term. The rest of the stations were in the range of 0.04 kg m −2 to 0.27 kg m −2 , which means that for most of them all GNSS solutions differ from each other by less than 3%. By comparing these results to the annual amplitudes from the RS, it can be noticed that various stations were characterized by various similarities to the RS. Firstly, there were stations for which the estimated values of the annual component on the basis on GNSS and RS were practically the same. DELF and 06260 may serve as an example. Obtained in this case the differences were usually below the value of the estimated for ERA-Interim accuracy (Table 4)

IWV Amplitudes
Because they took place in the mean IWV, the values of the annual oscillations did not significantly vary between the particular GNSS solutions (Figure 3) (Table A4 in Appendix A). More precisely, the biggest discrepancies were found for the MOPI station and amounted to 0.33 kg m −2 (GAM_DD_STP and GAM_DD_EL20) ( Table A5 in Appendix A). This value makes up about 4% of the total value of the annual component for this station and is the greatest deviation from all solutions both in terms of absolute and relative term. The rest of the stations were in the range of 0.04 kg m −2 to 0.27 kg m −2 , which means that for most of them all GNSS solutions differ from each other by less than 3%. By comparing these results to the annual amplitudes from the RS, it can be noticed that various stations were characterized by various similarities to the RS. Firstly, there were stations for which the estimated values of the annual component on the basis on GNSS and RS were practically the same. DELF and 06260 may serve as an example. Obtained in this case the differences were usually below the value of the estimated for ERA-Interim accuracy (Table 4), which was at the level of 0.02 kg m −2 . A similar situation occurred in the cases of BOR1, JOZE, WARE, and WTZR with corresponding values of 12425, 12374, 10410, and 10771. The annual amplitude values obtained from them mostly differ by less than the model accuracy and less than the annual amplitude estimation error, which was at the level of 0.06 to 0.09 kg m −2 . On the other hand, there were stations like EBRE and 08160, which annual amplitudes differ about 2.79 kg m −2 . For this pair of stations, we found also the biggest differences in the mean IWV value. However, next to EBRE, also GOPE, GRAZ, MATE, MOPI, and ZIMM stations where characterized by higher discrepancies w.r.t. RS stations. Similar to EBRE, each of them is at a long distance from the corresponding RS station, and, probably, this is the reason of the various annual signals. The only exception to this is the GRAZ station, which is located only 144 m higher than the corresponding 11240 station. The height difference between these two stations seems not to be big compared to the other mentioned cases. For example, the height difference between MOPI and 11035 is 335 m, and between MATE and 16320 it is 480 m. However, GRAZ and 11240 are situated in the Alps. Therefore, their local weather conditions are mountainous, which could cause the occurrence of the obtained differences despite the relatively small difference in height. The EBRE and 08160 stations varied from the mentioned cases, most probably due to the fact that, apart from the differences in height, they are also separated by big distance, so local weather conditions may also be completely different. The fact that the obtained results between GNSS and RS seem to be logical and explainable confirms the usefulness of GNSS in the investigation of changing prevailing weather conditions. Remote Sens. 2018, 10, x FOR PEER REVIEW 12 of 25 differences in the mean IWV value. However, next to EBRE, also GOPE, GRAZ, MATE, MOPI, and ZIMM stations where characterized by higher discrepancies w.r.t. RS stations. Similar to EBRE, each of them is at a long distance from the corresponding RS station, and, probably, this is the reason of the various annual signals. The only exception to this is the GRAZ station, which is located only 144 m higher than the corresponding 11240 station. The height difference between these two stations seems not to be big compared to the other mentioned cases. For example, the height difference between MOPI and 11035 is 335 m, and between MATE and 16320 it is 480 m. However, GRAZ and 11240 are situated in the Alps. Therefore, their local weather conditions are mountainous, which could cause the occurrence of the obtained differences despite the relatively small difference in height. The EBRE and 08160 stations varied from the mentioned cases, most probably due to the fact that, apart from the differences in height, they are also separated by big distance, so local weather conditions may also be completely different. The fact that the obtained results between GNSS and RS seem to be logical and explainable confirms the usefulness of GNSS in the investigation of changing prevailing weather conditions. Next to the annual amplitudes, we also estimated the values of semiannual signal for each of the analyzed solutions (Table A6 in Appendix). From all of the strategies, only those based on the PPP mode gave in the case of the SFER station a slightly higher value of semiannual amplitude than the rest of the analyzed time series. By comparing them to the RS (Table A7 in Appendix), it can be noticed that in the vast majority of cases differences between GNSS and RS were below 0.10 kg m −2 . Only in the case of EBRE and the mentioned SFER station, they were higher and amounted to 0.6 kg m −2 and 0.3 kg m −2 , respectively. Next to the annual amplitudes, we also estimated the values of semiannual signal for each of the analyzed solutions (Table A6 in Appendix A). From all of the strategies, only those based on the PPP mode gave in the case of the SFER station a slightly higher value of semiannual amplitude than the rest of the analyzed time series. By comparing them to the RS (Table A7 in Appendix A), it can be noticed that in the vast majority of cases differences between GNSS and RS were below 0.10 kg m −2 . Only in the case of EBRE and the mentioned SFER station, they were higher and amounted to 0.6 kg m −2 and 0.3 kg m −2 , respectively.

IWV Linear Trend
An analysis of the obtained results showed that there was a much smaller consistency in the linear trend value between the GNSS solutions (Table A8 in the Appendix A and Figure 4). Depending on the station, the obtained differences in linear trend varied from 0.015 kg m −2 /year, like in the case of BOR1, to 0.121 kg m −2 /year, like in the case of SFER. By complementing these values with the value of the linear trend estimation error (0.006-0.008 kg m −2 /year), it can be noticed that the mentioned differences were significant and not random. Therefore, they have to result (at least partly) from the adopted GNSS observation processing strategy. Similar to in the previous analysis, we compared the obtained results to the RS results. As it can be seen in Table A9 and Figure 5, the most convergent to the RS linear trend values were those obtained from the PPP solutions (BSW_PPP_VMF and BSW_PPP_GMF). Both in the case of applying GMF and VMF1, the absolute mean value of the differences between these two solutions and the RS solutions were at the level of 0.008 kg m −2 /year. Also, the selected values of the linear trend in the PPP solutions did not vary much from the RS trends. As an example, for the 12 stations in BSW_PPP_VMF, the linear trend value differs from the RS linear trend by less than its estimation error, whereas in the PPP_GMF such a situation occurred in the case of even the 14 stations. High similarity between both solutions and RS is also expressed by the lowest standard deviations (SD) of differences, which were at the level of 0.011 kg m −2 /year. For comparison to the PPP solutions, two other BSW solutions conducted in DD mode (BSW_DD_VMF and BSW_DD_GMF) had almost twice as high value of absolute mean difference, which amounted to 0.017 kg m −2 /year for the VMF and 0.016 kg m −2 /year for GMF solutions. The standard deviations were also two times higher than in the PPP solutions and amounted to 0.022 kg m −2 /year. By comparing the PPP solutions to the remaining differential solutions (all GAMIT), it can be also noticed that PPP presents more reliable results. The comparison of only the differential solutions looks interesting as well. Firstly, two solutions (BSW_DD_VMF and GAM_DD_VMF) conducted according to EPN AC (Analysis Centre) guidelines seems to be characterized by similar accuracy compared to the RS solutions. Despite several cases for which they had clearly different values of linear trend (e.g., HERS, SFER), they had the same absolute mean value of differences w.r.t. RS, which was 0.017 kg m −2 /year. SDs for these both solutions were not identical but were also very similar to each other and amounted to 0.022 kg m −2 /year and 0.023 kg m −2 /year for BSW_DD_VMF and GAM_DD_VMF, respectively. Secondly, applying the GMF/GPT combination instead of VMF1/ECMWF resulted in a slight improvement in the accuracy of the estimated linear trends. The absolute mean differences of the solution obtained on the basis of GMF/GPT were at the level of 0.016 kg m −2 /year. Thirdly, applying the GPT2 instead of GMF/GPT or VMF1/ECMWF seems to bring the best results of all differential solutions. Although GAM_DD_GPT2 still represented poorer accuracy than the PPP solutions, it is clearly visible that linear trend values from this solution were more consistent with the RS linear trends than those that took place in the case of GAM_DD_VMF1. As an example, the absolute mean differences from GAM_DD_GPT2 were at the level of 0.012 kg m −2 /year, whereas SD of these differences was at the level of 0.015 kg m −2 /year. This solution still contained some more significant discrepancies compared to the RS trends, but the highest one amounted to −0.035 kg m −2 /year (SFER). At the same time, in the BSW_DD_VMF and BSW_DD_GMF it was −0.063 kg m −2 /year and −0.061 kg m −2 /year (also SFER), whereas in case of GAM_DD_VMF it was 0.076 kg m −2 /year (REYK). Next to the using of GPT2, also, applying higher elevation mask brings benefits to the GNSS IWV solutions. As it can be noticed from the comparison of the GAM_DD_VMF and GAM_DD_EL20, there is a reduction of absolute mean differences value (to 0.014 kg m −2 /year) and SD of these differences (to 0.016 kg m −2 /year). Finally, as it was expected, the solutions obtained on the basis of the NIELL/STP combination (GAM_DD_STP) were characterized by poorer accuracy. The absolute mean difference between the linear trends from this solution and from RS was the lowest from all of the GNSS solutions and amounted to 0.022 kg m −2 /year.
A detailed analysis of the results showed also some aspects worth emphasizing like, e.g., differences in the case of SFER station for which BSW_DD_GMF and BSW_DD_VMF were characterized by much poorer accuracy (−0.063 kg m −2 /year) than GAM_DD_VMF (0.002 kg m −2 /year). Such clear discrepancies between two analogical solutions (BSW_DD_VMF and GAM_DD_VMF) occurred also in the cases of the REYK, JOZE, and HERS stations, but in these cases the solution from Bernese was more reliable. From the mentioned stations, HERS deserves special attention because it is only 4 km away from the 03882 station and the height difference between them is only 21 m. Due to the fact that it can be assumed that over both these stations there occurred similar weather conditions, the results obtained for GAM_DD_VMF are interesting. The opposite situation (higher accuracy of GAM_DD_VMF) can be found for MEDI, which is also located in close neighborhood of the corresponding RS station. Bernese was more reliable. From the mentioned stations, HERS deserves special attention because it is only 4 km away from the 03882 station and the height difference between them is only 21 m. Due to the fact that it can be assumed that over both these stations there occurred similar weather conditions, the results obtained for GAM_DD_VMF are interesting. The opposite situation (higher accuracy of GAM_DD_VMF) can be found for MEDI, which is also located in close neighborhood of the corresponding RS station. Due to the fact that both PPP solutions were deprived of such large and hardly explainable discrepancies, and were characterized by the highest consistency with RS, we decided to use them for a detailed comparison with RS. Figure 5 represents the linear trend value for BSW_PPP_VMF and RS. As can be clearly noticed from it, all the trends are positive, so in terms of trend character there is a 100% consistency between the GNSS PPP and RS method. This occurred even for stations for which we noticed the highest discrepancies in term of mean IWV or annual amplitude value, like EBRE with 08160 or MOPI with 11035. Generally, 65% of the stations were characterized by smaller discrepancies between GNSS_PPP and RS than the linear trend estimation error and 55% by less than the estimated mean ERA-Interim accuracy w.r.t. RS (from which we sourced meteorological parameters to the GNSS IWV conversion). The highest discrepancies between GNSS PPP IWV and RS IWV occurred for VILL and 08221. For the first of them, the linear trend value reached up to 0.016 kg m −2 /year, whereas for the second of them it was 0.044 kg m −2 /year. These stations are only 32 km away, and the height difference between them is 38 m. The estimated ERA-Interim accuracy for the 08221 station was at the level of 0.003 kg m −2 /year (Table 4). Consequently, nine times greater differences in the linear trend value cannot be explained by model accuracy. They also cannot be explained by the linear trend estimation uncertainty, which was at the level of 0.006-0.008 [kg m −2 /year]. Such big changes were probably not caused by real weather conditions, as well. A similar Due to the fact that both PPP solutions were deprived of such large and hardly explainable discrepancies, and were characterized by the highest consistency with RS, we decided to use them for a detailed comparison with RS. Figure 5 represents the linear trend value for BSW_PPP_VMF and RS. As can be clearly noticed from it, all the trends are positive, so in terms of trend character there is a 100% consistency between the GNSS PPP and RS method. This occurred even for stations for which we noticed the highest discrepancies in term of mean IWV or annual amplitude value, like EBRE with 08160 or MOPI with 11035. Generally, 65% of the stations were characterized by smaller discrepancies between GNSS_PPP and RS than the linear trend estimation error and 55% by less than the estimated mean ERA-Interim accuracy w.r.t. RS (from which we sourced meteorological parameters to the GNSS IWV conversion). The highest discrepancies between GNSS PPP IWV and RS IWV occurred for VILL and 08221. For the first of them, the linear trend value reached up to 0.016 kg m −2 /year, whereas for the second of them it was 0.044 kg m −2 /year. These stations are only 32 km away, and the height difference between them is 38 m. The estimated ERA-Interim accuracy for the 08221 station was at the level of 0.003 kg m −2 /year (Table 4). Consequently, nine times greater differences in the linear trend value cannot be explained by model accuracy. They also cannot be explained by the linear trend estimation uncertainty, which was at the level of 0.006-0.008 [kg m −2 /year]. Such big changes were probably not caused by real weather conditions, as well. A similar situation occurred in the case of GOPE with 11520 and GRAZ with 11240, for which no noticeable differences were found. It is difficult to assess which methods gave more reliable results; however, it is not insignificant that all these stations are at high altitudes. Perhaps there are some altitude-related aspects that caused this inconsistency.
Excluding these clearly unreliable stations, it can be noticed that rest of the differences did not have to result from numerical or environmental errors but rather represent some climatological information. As an example, linear trends obtained from the DELF and WARE stations are smaller than those corresponding to the 06260 and 10410 stations. However, if we take into account the fact that stations HERS and 03882 located to the west of them had an almost unnoticeable trend, and stations WTZR and 10771 located to the east of them had a higher trend, it can be assumed that the visible differences reflect real changes.
Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 25 situation occurred in the case of GOPE with 11520 and GRAZ with 11240, for which no noticeable differences were found. It is difficult to assess which methods gave more reliable results; however, it is not insignificant that all these stations are at high altitudes. Perhaps there are some altitude-related aspects that caused this inconsistency. Excluding these clearly unreliable stations, it can be noticed that rest of the differences did not have to result from numerical or environmental errors but rather represent some climatological information. As an example, linear trends obtained from the DELF and WARE stations are smaller than those corresponding to the 06260 and 10410 stations. However, if we take into account the fact that stations HERS and 03882 located to the west of them had an almost unnoticeable trend, and stations WTZR and 10771 located to the east of them had a higher trend, it can be assumed that the visible differences reflect real changes.

Discussion
The obtained results showed that there are differences in the analyzed GNSS IWV solutions. Generally, the short term variations of the GNSS IWV time series, like the annual amplitudes, were very consistent. The obtained differences w.r.t. RS were mostly caused by the area topography and not by the distance between the stations itself. As an example of this may serve EBRE with 08160 and WARE with 10410. In the case of the first pair of stations, which are located 147 km away from each other, we obtained the biggest differences from all analyzed stations, both in terms of annual amplitude and mean IWV value. On the other hand, stations WARE and 10410 are at 145 km distance and the differences obtained for them were much smaller. The reason of this inconsistency is the altitude difference between the EBRE and 08160 stations (Δh = 200 m) compared to WARE and 10410 (11 m). Despite this, it should be emphasized that the occurrence of the differences in seasonal components does not have to mean that there will be also discrepancies in terms of long-time changes. For example, the mentioned EBRE and 08160 stations for which the similarity between all GNSS solutions and the RS solution were at one of the highest levels.
All linear trends obtained from the analyzed solutions were positive, exactly like in the case of the estimated RS linear trends. These results were also compatible with other studies related to measuring long-term changes of IWV above Europe (for example, studies conducted by Bock et al.

Discussion
The obtained results showed that there are differences in the analyzed GNSS IWV solutions. Generally, the short term variations of the GNSS IWV time series, like the annual amplitudes, were very consistent. The obtained differences w.r.t. RS were mostly caused by the area topography and not by the distance between the stations itself. As an example of this may serve EBRE with 08160 and WARE with 10410. In the case of the first pair of stations, which are located 147 km away from each other, we obtained the biggest differences from all analyzed stations, both in terms of annual amplitude and mean IWV value. On the other hand, stations WARE and 10410 are at 145 km distance and the differences obtained for them were much smaller. The reason of this inconsistency is the altitude difference between the EBRE and 08160 stations (∆h = 200 m) compared to WARE and 10410 (11 m). Despite this, it should be emphasized that the occurrence of the differences in seasonal components does not have to mean that there will be also discrepancies in terms of long-time changes. For example, the mentioned EBRE and 08160 stations for which the similarity between all GNSS solutions and the RS solution were at one of the highest levels.
All linear trends obtained from the analyzed solutions were positive, exactly like in the case of the estimated RS linear trends. These results were also compatible with other studies related to measuring long-term changes of IWV above Europe (for example, studies conducted by Bock et al. (2014) [51], who on the basis of several techniques have also found only positive trends in Europe, or by Ning et al. (2016) [13], who for the analogous stations that we analyzed have found positive trends as well). Generally, the highest consistency between GNSS IWV and RS IWV long-term changes was found for PPP solutions. When both VMF1 and GMF were using in the PPP mode, the absolute mean difference of the linear trends w.r.t. RS and their standard deviations was the lowest of all of the analyzed strategies. In fact, most of the differences between the PPP solutions and the RS solution were below the linear trend estimation error. It should be kept in mind here that the first several years of the analyzed time span were characterized by poor satellites orbits and clocks quality, which are crucial in the PPP method, and yet the values of the linear trends seem to be more accurate than in the differential approach. This may indicate that the long-term troposphere changes obtained from GNSS are sensitive to network effects, which are unavoidable in the double difference approach.
With the exception of the PPP solutions, the best GNSS IWV solution was the one obtained on the basis of the GPT2 model (GAM_DD_GPT2). Its advantage over the solution in which GMF/GPT model (BSW_DD_GMF) was used results probably from the differences in the NWM (Numerical Weather Model) data that were used for mapping function coefficients estimation. In the case of GPT, it was monthly mean profiles from the ERA-40 from 1999-2002 period of time (23 pressure levels), whereas in case of GPT2 model, it was monthly mean profiles from ERA-Interim from 2001-2010 period of time (37 pressure levels). In addition, GPT2 is characterized by higher horizontal resolution and estimation of semi-annual terms of temporal variability (in the GPT it is only the mean and annual terms), as well as the phase (fixed in GPT). As it was confirmed by Lagler et al. (2013) [43], these advantages of the GPT2 model lead to the improvement of estimation of annual and semi-annual stations' height components from the radio-space geodetic techniques, for which proper estimation is correlated with proper estimation of the troposphere parameters [52]. On the other hand, it is hard to explain the advantage of this solution over the solutions that were based on the VMF1, since the empirical model like GPT2 does not take into account rapid and large weather changes. ZHDs sourced from ECMWF have been demonstrated to improve ZTDs and station height repeatability compared to the empirical models like GMF/GPT or NIELL/STP (Fund et al. 2010) [53]. However, these studies were conducted based on only one year of data from EPN, and although they confirmed the superiority of the GMF/GPT model over NIELL/STP, which was also reflected in our results, it is difficult to refer them to the investigation of long-time changes. Perhaps some effects related to the sensibility of VMF1 to the short-time weather changes are translated from one station to the other through baselines between stations. It would explain the lack of major differences between GMF and VMF1 in the PPP solutions and more visible differences when, instead of a weather prediction model, we used averaged and smoothed vales from ten years of data (a comparison of GAM_DD_VMF and GAM_DD_GPT2). Next to the mapping functions and ZHD a priori value, the elevation mask also plays an important role in the final GNSS IWV accuracy. Ning et al. [28] firstly demonstrated that its higher value has a positive impact on estimated long-term changes. Also, on the basis of our set of data we can confirm that applying a 20 • elevation mask instead of the routinely used 5 • increases the consistency between GNSS IWV and RS IWV. In our set of data this increase of quality was found in 75% of cases (comparison of GAM_DD_VMF and GAM_DD_EL20).
Finally, it has to be highlighted that beyond the imperfection of such indirect IWV measurements as used by us the GNSS IWV conversion formula (detailed discussed in [13]), obtained differences in linear trend values did not have to result from the differences in the GNSS and RS stations localizations or from the adopted GNSS processing strategy, but also from the inhomogeneity of the RS observations [54]. On the other hand, we do not have any other source of direct IWV observational data that would provide data for such a long time. Therefore, it seems reasonable to use the RS observations for the validation of indirect methods, but with the knowledge that some differences between them may result from a lack of RS metadata. However, the high consistency between RS IWV and GNSS IWV, especially in the case of the PPP approach, does not seem to be accidental, thus confirming the usefulness of GNSS measurements in climate studies.

Conclusions
In this paper we have analyzed 20 years of GNSS data that were reprocessed according to eight different processing strategies. Firstly, we have converted them to the IWV using the meteorological data sourced from the ERA-Interim. Secondly, we have validated them by using the IWV time series from neighboring RS stations. The comparative analysis concerned the mean IWV value, as well as annual and semi-annual amplitude, but mostly was focused on the linear trend. After analyzing the obtained results, we can state that: 1.
The highest consistency between the GNSS IWV linear trends and RS IWV linear trends was found for the PPP solutions. Both in terms of using VMF1 or GMF in PPP mode, estimated long term IWV changes were usually characterized by the highest accuracy.

2.
The highest differences between the GNSS PPP and RS were found only for these stations, which are located at higher altitudes. Therefore, some altitude-dependent error both for the GNSS PPP and RS methods should be considered in future works. 3.
DD solutions, widely regarded as the most accurate, were characterized by lower consistency w.r.t. RS than the PPP ones. The PPP solutions are realized individually, whereas the DD solutions are realized in regional network. This may indicate that, due to some network effects, the DD method may introduce to the troposphere solutions errors that affect the proper investigation of long term changes. 4.
GPT2 brings benefits to climate-related studies. This is probably due to the fact that this empirical model better represents long term tropospheric variations than numerical weather models like ECMWF, which are mostly focused on short-time forecasts.
Author Contributions: Z.B. and G.N. designed the research. Z.B. was responsible for long-term parameters estimation, results analysis, and writing the paper. G.N. was responsible for processing soundings observations and ERA-Interim data; he also took a part in the results analysis. M.F. was responsible for GNSS observations processing using Bernese GNSS Software. A.A. was responsible for GNSS observations processing using GAMIT software. All authors participated in the project and approved the final manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.  Table A3. Differences of mean value of IWV estimated from GNSS observations using different strategies (see Table A2) and sounding observations (see Table A1).   Table A5. Differences of annual amplitude of IWV (∆A A IWV ) estimated from GNSS observations using different strategies (see Table A4) and sounding observations (see Table A1).   Table A7. Differences of semiannual amplitude of IWV (∆A S IWV ) estimated from GNSS observations using different strategies (see Table A6) and sounding observations (see Table A1).   Table A9. Differences of linear trend of IWV (∆v IWV ) estimated from GNSS observations using different strategies (see Table A8) and sounding observations (see Table A1). Lowest differences marked by green, highest by red.