Comparison of GPS tropospheric delays derived from two consecutive EPN reprocessing campaigns from the point of view of climate monitoring

The main purpose of this research was to acquire information about consistency of ZTD (zenith total delay) linear trends and seasonal components between two consecutive GPS reprocessing campaigns. The analysis concerned two sets of the ZTD time series which were estimated during EUREF (Reference Frame Sub-Commission for Europe) EPN (Permanent Network) reprocessing campaigns according to 2008 and 2015 MUT AC (Military University of Technology Analysis Centre) scenarios. Firstly, Lomb–Scargle periodograms were generated for 57 EPN stations to obtain a characterisation of oscillations occurring in the ZTD time series. Then, the values of seasonal components and linear trends were estimated using the LSE (least squares estimation) approach. The Mann–Kendall trend test was also carried out to verify the presence of linear long-term ZTD changes. Finally, differences in seasonal signals and linear trends between these two data sets were investigated. All these analyses were conducted for the ZTD time series of two lengths: a shortened 16-year series and a full 18-year one. In the case of spectral analysis, amplitudes of the annual and semi-annual periods were almost exactly the same for both reprocessing campaigns. Exceptions were found for only a few stations and they did not exceed 1 mm. The estimated trends were also similar. However, for the reprocessing performed in 2008, the trends values were usually higher. In general, shortening of the analysed time period by 2 years resulted in a decrease of the linear trends values of about 0.07 mm yr−1. This was confirmed by analyses based on two data sets.

Abstract. The main purpose of this research was to acquire information about consistency of ZTD (zenith total delay) linear trends and seasonal components between two consecutive GPS reprocessing campaigns. The analysis concerned two sets of the ZTD time series which were estimated during EUREF (Reference Frame Sub-Commission for Europe) EPN (Permanent Network) reprocessing campaigns according to 2008 and 2015 MUT AC (Military University of Technology Analysis Centre) scenarios. Firstly, Lomb-Scargle periodograms were generated for 57 EPN stations to obtain a characterisation of oscillations occurring in the ZTD time series. Then, the values of seasonal components and linear trends were estimated using the LSE (least squares estimation) approach. The Mann-Kendall trend test was also carried out to verify the presence of linear long-term ZTD changes. Finally, differences in seasonal signals and linear trends between these two data sets were investigated. All these analyses were conducted for the ZTD time series of two lengths: a shortened 16-year series and a full 18-year one. In the case of spectral analysis, amplitudes of the annual and semi-annual periods were almost exactly the same for both reprocessing campaigns. Exceptions were found for only a few stations and they did not exceed 1 mm. The estimated trends were also similar. However, for the reprocessing performed in 2008, the trends values were usually higher. In general, shortening of the analysed time period by 2 years resulted in a decrease of the linear trends values of about 0.07 mm yr −1 . This was confirmed by analyses based on two data sets.

Introduction
Climate plays a key role in shaping the environment in which we live. It is a changing set of interconnected phenomena and, therefore, requires continuous research aimed at evaluating the current state of the atmosphere and predicting its future changes. Water vapour is one of the most important natural greenhouse gases. It is responsible for the Earth energy balance (Soden and Held, 2006) and it is one of the major factors that affects climate change. The exact extent of its impact on the heat-trapping quantity is still under discussion; however, it was already demonstrated that it is responsible for about 60-70 % of the Earth's surface temperature increase (Kiehl and Trenberth, 1997). Water vapour also plays a major role in shaping the dynamic processes in the atmosphere and the hydrological cycle. All these factors motivate scientists to monitor the variability of its content in the atmosphere. The amount of water vapour in the troposphere (integrated water vapour, IWV) can be measured by means of several techniques using radiometers, radiosondes, sun photometers and satellite sounders (e.g. AIRS, GOME), see e.g. Van Malderen et al. (2014). However, due to increased need for more climatological data, geodetic techniques (like global navigation satellite systems, GNSS; Doppler Orbitography and Radiopositioning Integrated by Satellite, DORIS; verylong-baseline interferometry, VLBI) are used for this purpose as well. Among them, GPS (global positioning system) observations started to play a significant role thanks to the availability of data covering 18 years. Moreover, GPS stations continuously provide data of high temporal and spatial resolution, in all weather conditions. In practice, the GPS signal propagation delay, caused by the physical properties of the neutral part of the atmosphere, is measured and exploited. Advanced processing of the GPS observations enables us to determine this delay with very high accuracy, below 2 mm (e.g. Bock et al., 2016).
The delay is commonly transformed to the zenith direction and is, therefore, called the zenith total delay (ZTD). It includes delays caused by the hydrostatic (ZHD, zenith hydrostatic delay) and the wet (ZWD, zenith wet delay) parts of the atmosphere. In principle, ZTD reflects the state of the troposphere, so it can be used together with selected meteorological parameters (surface pressure and temperature) to estimate IWV (Bevis et al., 1992). Then, using long-term observations, it is also possible to study its changes over time (e.g. Emardson et al., 1998;Nilsson and Elgered, 2008;Wang and Zhang, 2009), which constitutes the beginning of its use in climate research. However, it is worth noticing here that IWV obtained from ZTD and meteorological parameters is affected by some errors. They are a result of the uncertainties in (i) the ZTD estimation, (ii) the ZHD modelling, (iii) the method of interpolation (in all those cases when meteorological parameters have to be interpolated e.g. from weather models) and (iv) the ZWD to IWV conversion formula (e.g. Hagemann et al., 1992;Bock et al., 2007;Van Malderen et al., 2014;Ning et al., 2016). Therefore, several studies have been conducted to determine the best methods of reducing these errors and uncertainties. These methods focused mainly on the homogenisation of the ZTD time series (e.g. Hagemann et al., 1992;Bock et al., 2014) or using various methods of interpolation (Wilgan et al., 2014). It is worth keeping in mind that GPS achieved full operational capability in 1995. The processing strategy has been changed several times since then. As in other climate studies, the analysed data should also be as homogeneous as possible (Bengtsson et al., 2004). In this context, using reprocessed data will minimise errors resulting from a switch in the processing strategy. Nowadays, their potential as a source of information for climate research is constantly growing (e.g. Pacione et al., 2014). Moreover, the adopted processing strategy could also affect the estimated ZTD values. This is particularly important when changes in the water vapour content are very small and any mismodelling of the tropospheric parameters could also distort the estimated trend. As already presented by Nilsson and Elgered (2008) or Baldysz et al. (2015), extending the analysed period even for 2 additional years can change the nature of the trend. This confirms the need for long time series data for climate research. ZTD shows correlation with time-dependent temperature (Guerova, 2013) or water vapour content (Yong et al., 2008). Its size and seasonal variability are related to such factors such as latitude and altitude or distance to large masses of water. Therefore, it could also provide information about the prevailing local weather conditions (Jin et al., 2007). Having these facts and taking into account that ZTD is a direct GPS product, we decided to focus on this parameter. To investigate the impact of an adopted GPS processing strategy on ZTD parameter and, therefore, on its climatological applications, we compared two different data sets. The first of them consisted of the results of the MUT (Military University of Technology) tests conducted before the first EPN (EUREF Permanent Network;Bruyninx, 2004) reprocessing project, extended to 2013 according to the same strategy. MUT results are called Repro1; they were obtained on the basis of full EPN reprocessing made by MUT and did not play a role in the official EPN Repro1 project solution Kenyeres et al., 2009. These data were already analysed by Baldysz et al. (2015). The second one was the MUT's official contribution to the latest EPN reprocessing campaign (called here Repro2). As the ZTD parameter should reflect the real amount of delay caused by the troposphere, theoretically estimated ZTD during both reprocessing campaigns should be the same. Consequently, the designated oscillations and linear trends should also be identical for time series from exactly the same period of time (from the Repro1 and Repro2 campaigns). However, the use of various software and computing strategies, including models and parameters, may influence the values. For this reason we decided to analyse the differences between these two solutions. Analyses were conducted for the ZTD time series of two lengths: a shortened 16-year period to maximise the spatial coverage (57 stations), and a full 18-year period to ensure the maximum length of the analysed period of time (28 stations).
This paper is structured as follows. Sections 2 and 3 describe the data and methods used for analysis. Sections 4 and 5 are dedicated to the results of short-term (seasonal components) and long-term (linear trends) analyses. The obtained results are discussed in a broader perspective in Sect. 6. The last section provides the summary and the conclusions of the performed analysis.

Analysed data
The ZTD is a consequence of electromagnetic wave delay (T ) which is caused by refraction and attenuation in the troposphere. This delay T is defined by the following formula (Bevis et al., 1992): where c is the speed of light in vacuum, τ is the delay measured in the unit of time and N is the neutral atmospheric refractivity (Davis et al., 1985): where k i (i = 1, 2, 3) are constants, R d is a specific gas constant for dry air, ρ d is the total mass density of the atmosphere, P w is the partial pressure of water vapour, Z w is the compressibility factor and T is the temperature in Kelvin.
Atmos. Meas. Tech., 9, 4861-4877, 2016 www.atmos-meas-tech.net/9/4861/2016/ The delay caused by the hydrostatic part of the atmosphere is given by the integral of the first term of Eq. (2) while the influence of the wet part of the atmosphere is given by the integral of the remaining two terms of Eq. (2). Most GPS satellites which are tracked by receivers for observations are not in the zenith direction. Therefore, the value of delay between each pair of a satellite and an antenna is calculated not in the zenith but in a slant direction (STD, slant tropospheric delay). Hence, it depends on the satellite's zenith distance. For relating STD to ZTD, it is necessary to use mapping functions which are approximately equal to 1/ sin e, where e is the elevation angle. However, for precise position determination, this approximation is not sufficient and the following function given by Marini (1972) and normalised by Herring et al. (1992) has to be used: where e is the elevation angle and a, b, and c are coefficients which are related with the state of the atmosphere (determined differently for different mapping functions). Mapping functions are used separately for hydrostatic and wet parts of delay. The total value of delay in the zenith direction can thus be expressed as follows: where mf(e) hyd is a mapping function for hydrostatic delay and mf(e) wet is a mapping function for wet delay, G n and G e are north and east troposphere gradients (e.g. MacMillan, 1995;Chen and Herring, 1997;Bar-Sever et al., 1998), and α is the azimuth. The ZHD in the above formula is usually adopted from a model, while ZWD is one of many unknowns which are estimated during position determination. For this reason, it may be affected by adopting various calculation strategies (e.g. Ning and Elgered, 2012). Therefore, obtaining the most reliable value of ZTD during the processing of GPS observations is one of the major tasks for scientists who focus on atmospheric research. In this paper the MUT's individual contributions to the EPN Repro1 test (non-official contribution covering the whole of Europe) and EPN Repro2 (official contribution covering the whole of Europe) projects were compared. This enables us to identify how updating the GPS processing strategy affects the tropospheric parameters. The EPN (www.epncb.oma.be) is a network of permanent GNSS stations built on the basis and as a densification of the IGS global network (International GNSS Service) in Europe. Together with the UELN (United European Levelling Network), it plays a role of the EUREF key infrastructure. The EUREF (www.euref.eu) is the IAG (International Association of Geodesy) sub-commission and its main goal is to define, realise, maintain and provide the European Reference Frame (currently ETRS89 and ETRF2000). The EUREF is a joint voluntary effort of many research agencies and national mapping and cartographic agencies. The EPN has been operated continuously since 1996. During this time, various algorithms, models, parameters and software were involved to process GPS data collected by the stations. This resulted in collecting and archiving inhomogeneous sets of solutions which, in consequence, disabled proper analyses of long time series. Such a state of affairs led to decisions in 2008 and 2013 concerning recalculations of all data (from the level of observation files) according to one coherent strategy reflecting the current capabilities. This task was done inter alia by the MUT, one of the EPN Analysis Centres (MUT AC) in the frame of the special EUREF project called "EPN reprocessing". The first calculations Figurski et al., 2009;Sohne et al., 2010) covered all EPN stations which were operated in the EPN network from January 1996 through to December 2007. The second campaign (2015) covered all EPN stations which were operated from January 1996 to December 2014. To align the analysed periods, we decided to recalculate the data from January 2008 to December 2013 according to the Repro1 strategy. This ensured maximum coherence for the comparison of the solutions from these two reprocessing campaigns.
The knowledge and capabilities of physical modelling changed in the years between the Repro1 and the Repro2 campaigns. Consequently, calculation strategies in both reprocessing campaigns were not identical and could result in differences in ZTD values. The first major difference between the MUT's contributions to Repro1 and Repro2 was in the applied software. Bernese 5.0 software (Dacht et al., 2007) was used in Repro1 while GAMIT 10.5 software  was used for the Repro2 calculations. Both solutions were based on the GPS system only and the network approach. Reprocessed IGS orbits (IGS repro1) were used in Repro1 while reprocessed orbits provided by CODE in 2013 were used in Repro2. This implies that the solutions were expressed in different reference frames (Repro1 in IGS05 and Repro2 in IGb08). Consequently, antenna models were also changed. Besides the alignment to the specific reference frame, the source of PCC (phase centre correction) was also different. IGS type mean and EPN individual calibrations were used in Repro1 while only the IGS type mean calibrations were used in the Repro2 campaign. In terms of geophysical phenomena modelling, only the ocean tidal loading model (based on FES2004, Finite Element Solution 2004) was applied to Repro1. In Repro2, the FES2004 as well as atmospheric tides based on the ECMWF CMT (European Centre for Medium Range Weather Forecasts convective momentum transport) and atmospheric non-tidal loadings based on the NCEP (National Centre for Environmental Prediction) were taken into account. These models were applied to data processing at the observation level, as described by Tregoning and Herring (2006). Along with geophysical phenomena, the impact of ionosphere on phase measurements was also removed in various ways. In the first reprocessing, 4864 Z. Baldysz et al.: Comparison of GPS tropospheric delays only the first-order ionospheric effects were eliminated by using a linear combination. In the second one, all three orders of ionospheric corrections were modelled according to Petrie et al. (2010). Another distinction lies in the adopted elevation mask: 3 • in Repro1 and 5 • in Repro2. Moreover, during the first reprocessing campaign ZTD solutions were fixed to weekly coordinates while during the second one they were fixed to daily coordinates. The differences mentioned above are not directly related to the troposphere modelling. However, the accuracy of GPS-derived ZTD is affected by elevation-angle-dependent errors which could be caused by atmospheric mapping functions (Stoew et al., 2007), antenna phase centre variations (Schmid et al., 2007) or signal multipath. (Elósegui et al., 1995). For factors directly related to the state of troposphere, the biggest differences between Re-pro1 and Repro2 occurred in the mapping functions and the a priori adopted ZHD. In Repro1, the Niell mapping function (NMF; Niell, 1996), was used with coefficients for hydrostatic and non-hydrostatic components obtained from radiosondes profiles. Another approach was adopted in the second reprocessing campaign. Instead of NMF, Vienna Mapping Function 1 (Boehm et al., 2006), was used with coefficients for hydrostatic and non-hydrostatic components delivered by the ECMWF (European Centre for Medium Range Weather Forecasts). The estimation of the influence of various mapping functions on troposphere modelling was the goal of several studies (e.g. Tesmer et al., 2007;Vey et al., 2006a). However, these studies were mostly focused on coordinates or short-term ZTD time series and, therefore, it is difficult to draw conclusions about their impact on climate applications. The adopted processing strategies for both reprocessing campaigns are summarised in Table 1.
The two sets of data described above were derived from 57 EPN stations which had been operated continuously since at least 1998. In the case of the EPN network, the longest possible time series come from the stations which were launched simultaneously at the beginning of 1996. However, there are only 30 such stations. This small network does not provide adequate spatial distribution which is necessary for analyses of temporal and spatial changes of tropospheric parameters in such a large area. Therefore, we decided to analyse 16-year ZTD time series for both reprocessing campaigns (January 1998 to December 2013), thereby doubling the number of involved stations in our analysis. The time series of stations already in operation before 1998 were limited to January 1998 to ensure maximum homogeneity. Full 18-year ZTD time series were also studied. This approach, similar to in our previous study (Baldysz et al., 2015), enables us to investigate whether the introduction of two additional years has the same impact on the data from both reprocessing campaigns.
Besides the length of the ZTD time series, their quality is particularly important for climate studies, as the investigated changes can be very small. The data quality should not be disturbed by gaps in time series and a large number of outliers. Therefore, ZTD screening was conducted in both solutions sets according to approach described by For all stations, the percentage of all rejected solutions was at the level of 0.42 %. The largest number of rejected solution occurred for the SFER station (Spain) and amounted to 4.18 %. Furthermore, a data availability analysis was performed -real numbers of data were compared to theoretical numbers of data (theoretical maximum number of hourly solutions available from a station). Stations with less than 90 % of the theoretical number of data were rejected. This requirement had to be met by each station in both sets of data (from the Re-pro1 and Repro2 campaigns). Only two of all the analysed stations, ANKR (Turkey) and SVTL (Russian Federation), had to be excluded from further analysis. The final step in preparing the data for further analysis was manually removing discontinuities by using site log files available for all stations. They were related to the antenna or receiver changes; therefore they occurred at the same time in both sets of data. We found discontinuities in every station. On average, 3.12 discontinuities were found in the ZTD series of station in our data set. The last segment of the time series was chosen as the reference to which the other segments were shifted. This approach was used in both data sets. The spatial distribution of stations included in the analysis is shown in Fig. 1.

Analysis of ZTD time series
Seasonal variations in the ZTD time series are important in climate analyses. On one hand, they carry information about prevailing weather conditions in a given region. On the other hand, assuming their occurrence is important when the linear trend is determined using the least squares estimation (LSE) approach. To ensure the homogeneity of the compared re-Atmos. Meas. Tech., 9, 4861-4877, 2016 www.atmos-meas-tech.net/9/4861/2016/ sults (from the Repro1 and Repro2 campaigns), the authors adopted the same methodology as in the case of our previous work (Baldysz et al., 2015) where data from the Re-pro1 campaign were analysed. The only difference lies in the temporal resolution of the used data: in our previous work hourly data were analysed, whereas in this work, we chose to work with daily data for both data sets. Daily values were applied instead of hourly ones, due to the fact that hourly data are more highly correlated; therefore they could disturb a proper trend analysis. Information concerning oscillations in the ZTD time series was obtained by preparing the Lomb-Scargle periodograms for every station. According to this method, the frequency spectrum is estimated by fitting the linear squares of sine and cosine models to the observed time series (Lomb, 1976): where x(t i ) is the observed time series at time t i , a and b are constant amplitudes, ω is the angular frequency, is the additional phase (required for the orthogonalisation of the sine and cosine model functions when the data are unevenly spaced) and n i is the noise at time t i . Based on this method, the main frequencies in the ZTD time series were found. The periodograms were prepared for both sets of data (Repro1 and Repro2). Similarly to our previous work, we removed the annual oscillation from the signal to investigate the smaller amplitude oscillations. Waveforms of periodograms obtained for the same stations but different reprocessings, match each other almost perfectly. Therefore, we present results for the most recent studies (Repro2). Figure 2a presents examples of periodograms for five stations: GRAS (France), GLSV (Russian Federation), SFER (Spain), MAS1 (Spain) and RAMO (Israel). The results of the seasonal component determination for the whole signal spectrum are presented in Fig. 2a (top). The periodograms for the same stations for the signal after removing the annual period are presented in Fig. 2a (bottom). Every station had a strong annual signal which is caused by the annual weather oscillation (the highest temperature in summer and the lowest temperature in winter) characteristic for midlatitudes. For the majority of stations, a distinct semi-annual oscillation was also noticed. For some of the analysed stations triannual (CASC, SFER, VILL, MAS1) and even quadrennial signals were found (SFER, MAS1). All these stations are located in the south-western part of the analysed network in Spain (CASC, SFER, VILL) or in the Canary Islands (MAS1), so probably the occurrence of some Atlantic currents affected them in such way. As well as these oscillations, characteristic oscillations with about a 640-day period and about 5.0 mm amplitude (Fig. 2b) were found for a group of stations located in northern Scandinavia (TRO1, KIRU, KIR0, SODA, VIL0). Their occurrence can be caused by various phenomena. However, we assume that they are probably related to a similar geographic location but this effect is still under investigation. The periodograms are presented for periods of not more than 2 years in order to improve their readability. They were originally prepared for investigating lower frequencies (even 4 years); however, such frequencies were not observed in the ZTD time series over Europe, or they were not sufficiently clear. The trend value was estimated by means of the LSE method taking into account the results obtained from the Lomb-Scargle periodograms (annual, semi-annual, triannual and quarter-annual oscillations were fitted together with the linear trend in one matrix) and white noise assumption. We also carried out first tests about the influence of using another noise model (autoregressive noise of the fourth order (AR(4)) on the estimated trend values and uncertainties during the linear trend detection. Based on our time series, we obtained linear trend values that were almost exactly the same for both assumptions. The investigated differences between these two noise models occurred for trend uncertainties, which were about 2 times higher for the AR(4) model. However, more important is that the variability of the ZTD linear trend uncertainties (between various stations) obtained on the basis of www.atmos-meas-tech.net/9/4861/2016/ Atmos. Meas. Tech., 9, 4861-4877, 2016 the AR(4) model do not differ much from the variability of the uncertainties obtained on the basis of only white noise assumption. At this step, it is hard to estimate how accurately the AR(4) model is able to describe signals related to the tropospheric changes, next to software noise or incorrectly modelled geophysical oscillations and, therefore, assess the nature of the resulting larger trend uncertainty values. Further application of the AR(4) model requires additional studies, e.g. analysis regarding the stationarity of the ZTD noise, which will be the goal of our future work. Therefore, at this step we decided to apply only white noise assumption for considerations presented in this paper. The values of the linear trends determined by the LSE approach were very small for some stations; therefore we performed the Mann-Kendall trend test (Mann, 1945;Kendall and Stuart, 1970) to test the significance of a trend in the time series. All calculations were conducted according to the following formula: where X i andX j are the time series while i = 1, 2, 3n − 1 and j = i + 1, i + 2, i + 3, .., n. S is the statistical factor for which the sign, positive or negative, describes the character of linear trend. The assessment of the trend line significance is tested by a statistical test in which S is transformed to the standard normal test statistic. We adopted a significance level α = 0.05 here.

Results of seasonal analysis
Seasonal analyses were performed for the 16-year ZTD time series to ensure higher density in spatial distribution. As mentioned in the previous section, both reprocessing campaigns gave very similar results. In more detail, the average value of the annual oscillation (considering all the stations shown in Fig. 1) was 46.8 mm from Repro1 and 46.4 mm from Repro2. In both data sets, the highest amplitude of the annual period was found for the same station (TORI, Italy). Its value was slightly different for both reprocessing campaigns (63.3 mm in Repro1 and 63.1 mm inRepro2). The lowest value of the annual oscillation was found for the RAMO station (Israel) and it was 13.8 and 14.1 mm for Re-pro1 and Repro2 respectively. Generally, the differences between the estimated amplitudes of the annual periods were lower than 1 mm, except for four stations: DELF (the Netherlands), HOFN (Iceland), TERS (the Netherlands) and WSRT (the Netherlands), where the differences amounted to 1.1, 1.2, 1.3 and 1.2 mm respectively. The highest amplitude of the semi-annual periods was also found for the same station (JOEN, Finland) in both campaigns. Its value differed slightly and was 11.9 mm in Repro1 and 12.1 mm in Re-pro2. There was also no significant difference in phase for the seasonal component. For each of the 52 stations, the maximum of ZTD occurred on the same day of the year for both campaigns. A small phase shift (not more than 4 days) was found for only five stations: DELF (the Netherlands), MAS1 (Spain), SODA (Finland), VILL (France), and VIL0 (Sweden). Even then, the mean value of the maximum ZTD was equal in both campaigns and occurred on the 214th day of the year. Due to the strong similarity of the results, only Repro2 is shown in Fig. 3. In Europe, the maximum temperature, which is correlated with water vapour, occurs in late July and early August. However, it can be clearly seen that there are differences between regions. In the Mediterranean area, the maximum value of ZTD occurred on a later date than in areas in northern Europe. This is probably due to the fact that during the summer the sun is higher above the horizon in the regions closer to the Tropic of Cancer. In these areas, the sun supplies energy for the longest time period (compared to the regions surrounding the Arctic Circle), intensively heating them. All results are presented in Table A1 in Appendix A. The information obtained from seasonal analyses, supplemented by the mean ZTD value, could be useful for investigations of prevailing weather conditions, because of the correlation between the temperature and the ZTD value. We give some examples here. The southern part of Europe, which is within the influence of the Atlantic Ocean and the Mediterranean Sea, is located in the humid subtropical climate zone. It is clearly visible that ZTD time series from stations like CASC (Portugal), SFER (Spain), MAS1 (Spain), or CAGL (Italy) are affected by prevailing weather conditions which are typical for the neighbourhood of warm water reservoirs. This is mainly reflected in the smallest annual oscillations in Europe (e.g. CASC 22.5 mm, SFER 24.4 mm) and the largest mean ZTD values (e.g. CASC 2.425 m, SFER 2.422 m). RAMO (Israel) is also located near the Mediterranean coast. However, despite the proximity of warm water masses, it is characterised by much lower humidity and annual amplitudes (the annual amplitude at the RAMO station reaches 13.8 mm and the mean ZTD value is up to 2.153 m). This can be ascribed to the location of the station, at an altitude of 893.1 m above sea level, and on the border of the moderate Mediterranean and dry arid climate zones. In the central part of Europe it is possible to distinguish less (in the interior of the continent) and more (closer to the ocean) humid vari-eties of the temperate climate zone. Stations located within the humid varieties usually have higher mean ZTD values and smaller annual amplitudes compared to the stations situated in the interior of the continent. This can be illustrated by examples of stations like HERS (UK) with 2.406 m of the ZTD mean and 41.0 mm of the annual amplitude or WARE (Belgium) with 2.374 m of the ZTD mean and 44.7 mm of the annual amplitude compared to the LAMA (Poland) with 2.359 m of the ZTD mean and 52.3 mm of the annual amplitude or PENC (Hungary) with 2.351 m of the ZTD mean and 52.8 mm of the annual amplitude. The different characters of the ZTD time series from a few selected stations are presented in Fig. 4. Of course, the topography near the station should also be taken into consideration when assigning features which are characteristic of selected climate zones to the ZTD time series. Mountain massifs might affect the formation of mountainous types of climate, because they disturb the movements of air masses, thereby forming e.g. precipitation shadows.

Results of trend analysis
The trend values for 57 EPN stations were estimated using the LSE method, and in the next step their statistical significance was confirmed using the Mann-Kendall trend test (for both reprocessing campaigns). Only stations with confirmed trend were used to determine the average trend for the whole of Europe. Trend values were estimated based on daily data and were in the range of −0.70 to 0.60 mm yr −1 .
The corresponding errors were in the range of 0.06 to 0.11 mm yr −1 . In Repro 1, the occurrence of a trend was confirmed for 53 of 57 stations. It was not confirmed for CASC (Portugal), MATE (Italy), MEDI (Italy) and TORI (Italy) stations. For the 53 stations the average value of the trend was equal to 0.12 mm yr −1 . 34 stations had positive trends and 19 stations had negative trends. In Repro2, the Mann-Kendall trend test found no trend for five stations: GRAS (France), GRAZ (Austria), MATE (Italy), SFER (Spain) and TRO1 (Norway). For the remaining stations, the average trend value was 0.04 mm yr −1 , 30 stations had positive and 22 stations had negative trend values. Results for all stations included in the analysis are presented in Fig. 5.
Stations without a significant trend according to the Mann-Kendall test are indicated with black dots. The geographical distribution of the trends is similar for both cases of processings. The strongest discrepancies occur in the values of the negative trends, which are noticeably larger in Repro 2 than in Repro1. This is also reflected in the lower average trend value for the whole of Europe for Repro2. The area of north-eastern and south-eastern Europe is characterised by positive trend values on both maps, with some discrepancies in the values. Similar consistency of the trends occurs in the area of western and south-western Germany, the Netherlands, Belgium and north-eastern France, where only   negative trends occur in both the Repro 1 and the Repro2 campaigns. A clear coherence of trends for such large areas proves that they were not a result of random phenomena, but they rather followed consistent changes. The strongest geographical inconsistency in the trend pattern takes place in central and southern Europe. In this area, trends with the opposite signs exist in close neighbourhood (for the Repro1 and Repro2 campaigns). However, it is worth noticing that there is a large topographical variety (the presence of highlands and mountain ranges) which has a very large influence on the movement of air masses and, therefore, on the formation of the prevailing weather conditions.
To obtain more precise information about the extent of the differences for each station, the trend value derived from the first reprocessing was subtracted from the trend value ob-tained from the later reprocessing. In Repro2, 42 stations had smaller trend values than in Repro1, 13 stations had higher values and only 2 stations had exactly the same value of the trend (KIR0, Sweden and SODA, Finland). The average value of the trend differences is 0.07 mm yr −1 with a maximum trend difference of 0.25 mm yr −1 . More details concerning the differences in the trend values between the Re-pro2 and the Repro1 campaigns are shown in Fig. 6. Now, we constrain our considerations to the 28 EPN stations, which have been operated continuously since 1996. In Repro1, the Mann-Kendall trend test gave no trend for 2 stations: DELF (Netherlands) and WARE (Belgium). For the rest of them, the mean value of the trend was 0.19 mm yr −1 , 5 stations had negative trends and 21 had positive trends.   Fig. 7. These maps clearly show why the mean value of trends for the 18-year ZTD time series is higher than for the 16-year ZTD time series. This is because most stations with negative trends (which occurred especially in the area of western and south-western Germany, the Netherlands, Belgium and north-eastern France) were not included in this analysis (too short time series). This resulted in greater consistency and homogeneity in the estimated ZTD trends because most stations which were included in the 18-year ZTD analysis had positive trends. For the 18-year time period, in Repro 1, only two stations had no trend identified by the Mann-Kendal test. In Repro2, the Mann-Kendall test found six stations without significant trend. Therefore, the spatial distribution of stations in Repro2 looks poorer than in Repro1. It is also worth noticing that for the 18-year period the Mann-Kendall trend test gave no significant trend for other stations than for the shorter period. This is particularly noticeable for  results based on the data from the Repro1 campaign. For this set of data and shorter time period, no significant trends were found for different stations compared to the longer period.

Atmos
In the Repro2 campaign, no significant trend was found for only two stations for both time series lengths, and they were GRAS (France) and MATE (Italy). However, in both cases in central and southern Europe, the negative trend areas occurred next to the positive ones, which is consistent with the results obtained from analysis of the 16-year ZTD time series. A similar consistency in the geographical distribution of trends between the shorter and longer periods of time can also be seen in the area of north-eastern Europe where only positive trends occurred. Finally, trends obtained from Re-pro2 were smaller than those from the Repro1 campaign, regardless of the length of the analysed time period, confirming the impact of the GPS processing strategy on the ZTD trends.
To obtain information about the extent of the differences in the trend values for each station, the trend value derived from Repro1 was subtracted from the trend value obtained from the second campaign. 25 stations had a lower value of the linear trend in Repro2 compared to Repro1 and only 3 stations had a higher value. None of the stations had exactly the same trend value. All results related to the analysis of the 18-year ZTD time series are given in Table 3 in Appendix. The details of the differences between both reprocessing campaigns are shown in Fig. 8.
Differences in the number of stations for which the trend was confirmed (for both lengths of the time series) and differences in the ZTD trend values confirm that it is necessary to work on exactly the same period of time to analyse the changes of the spatial distribution of the ZTD trends. Two additional years in the length of the time series may be important not only for the estimated value of the trend but also for its sign. In both campaigns there were differences in trend values and their signs between shorter and longer periods of time (e.g. WARE station in Repro1 and WTZR station in Re-Atmos. Meas. Tech., 9, 4861-4877, 2016 www.atmos-meas-tech.net/9/4861/2016/ pro2). However, both in Repro1 and Repro2, the values of the trends for the longer periods were generally larger than for the shorter ones. It should also be emphasised that the values of the differences between these two periods were at the same level and they had the same characteristics, both for Repro1 and Repro2, as shown in Fig. 9.

Discussion
The performed comparative analysis showed a clear differences in the obtained results. These differences were investigated both for the used data set (from various reprocessing campaigns) and the length of the time series (16-year time series and 18-year time series). They mainly concerned the values of the linear trends. The linear trend uncertainties were similar for the same time span, regardless of the adopted processing strategy. A slight decrease in their value was observed for the longer time period. Generally, the trend values in Repro2 were smaller (for both lengths of the time series) than in Repro1. These differences were significant for about half of the stations, for which discrepancies between two solutions were higher than ZTD trend estimation uncertainties. As in the case of the shorter time period, 28 out of 57 stations had differences in trend values between both reprocessing campaigns in the range of 0.09 to 0.25 mm yr −1 . For the longer time periods, these differences were significant for more than 60 % of stations and they also varied in a similar range (0.08 to 0.26 mm yr −1 ); however it should be kept in mind that, for this length of the time series, the corresponding trend uncertainties were smaller than for shorter ones. This may mean that differences arising from the use of different GPS processing strategies become more significant for longer periods of time. In view of the fact that in both reprocessing campaigns the same raw data were used, these differences should be ascribed to the processing strategy.
To assess the possible impact of the mapping function on the long-term time variability of the ZTD parameter, we selected 28 stations, for which we first analysed the ZHD time behaviour, which is responsible for about 90 % of the total value of the delay. The ZHD value was calculated from VMF1, because NMF, which was used in Repro1, is dependent only on the day of the year, latitude, and altitude above sea level. Therefore, the NMF does not have a time dependence over timescales of many years. The mean trend value of ZHD from VMF1, for all 28 selected stations, was −0.07 mm yr −1 , while the mean difference in ZTD trend value between Repro2 and Repro1 (for the same stations) was −0.09 mm yr −1 . Of course, at this stage it cannot be clearly stated that the adopted a priori ZHD value could have such a significant influence on the ZTD trend value. However, it is worth noticing that the negative trend value of the ZHD from VMF1 generally agrees with the smaller value of the ZTD trends obtained from the Repro2 campaign. There was only 0.02 mm yr −1 difference between the ZHD and the VMF1 mean trend value and the differences between Repro2 and Repro1. Nevertheless, this similarity does not completely solve the problem of the discrepancies in trends and we cannot conclusively state that the mapping function is responsible for the major part of the differences values. A detailed analysis for individual stations showed, in fact, that the trend value adopted a priori ZHD data were not equal to the differences: Repro2 -Repro1. For most stations, these two parameters have similar signs, and some of them even have nearly identical values. For some stations, such as DELF (Netherlands), MAS1 (Spain), METS (Finland), ONSA (Sweden), RIGA (Latvia) or VILL (Spain), the adopted a priori ZHD cannot explain the differences in trend values between the Repro2 and the Repro1 campaigns (see Fig. 10).
These two contributions of MUT into the reprocessing of EPN were different, not only in case of the applied tropospheric model. Therefore some other factors should be taken into account for finding out the reasons for the differences in the ZTD trends between both cases of reprocessing. Several studies were carried out in the field of the influence of GPS processing strategies on final GPS products. In most cases, they were focused both on the ZTD parameter and on vertical coordinates, due to the fact that they are correlated (Rothacher, 2002). Tregoning and Herring (2006) showed that a priori ZHD value had influence on both coordinates (vertical) and the ZTD value. Thomas et al. (2011) analysed, besides using various sources of meteorological data and mapping functions in troposphere modelling, the influence of the applied ATML (atmospheric loading model) and the type or individual antenna calibration on the ZTD parameter. They found out that the differences in ZTD value between these solutions reached up to 10 mm. The impact of the satellites cut-off angle was investigated inter alia by Vey et al. (2006b), who studied differences in the PW (precipitable water) obtained from GPS and found out that switching elevation mask from 5 to 15 • during GPS observation processing caused differences in PW value of about 0.75 mm. All these studies, however, were mostly focused on short time series (up to 1 year). Therefore their possible impacts on the secular variations of the ZTD parameter are difficult to estimate and there is a necessity to extend the studies to long time (even 20 years) series. An analysis of such a long-term dependence between elevation cut-off angle and trends in atmospheric water vapour content was carried out by Ning and Elgered (2012). They used 14 years of observations from 12 GPS and 7 radiosonde sites and processed them using eight different elevation masks (from 5 to 40 • ). The best agreement between IWV (integrated water vapour) from GPS and from radiosondes was obtained for 25 • elevation mask. Those studies point out that the mapping function has an influence on linear trends. However, differences in ZTD linear trends resulting from switching elevation mask from 3 to 5 • are still unclear and have to be explicitly defined due to the fact that these two cut-off angles are the most commonly used during official reprocessing of EPN.
www.atmos-meas-tech.net/9/4861/2016/ Atmos. Meas. Tech., 9, 4861-4877, 2016  Besides the differences of trend values between Repro2 and Repro1, it is worth noticing that both reprocessing campaigns gave almost the same results for differences between longer and the shorter periods of time (Fig. 11). This means that the adopted processing strategy did not directly influence the detected changes caused by the additional 2 years. Therefore, it could be indicated that the effect of non-linear changes in the state of the troposphere (caused by weather conditions which occurred during the additional time period) had the same impact on the trend value, regardless of the method of ZTD estimation.

Summary and conclusions
The ZTD time series obtained from two different EPN reprocessing campaigns, processed by the MUT, were analysed in this paper. The purpose of this work was to assess the differences in the time variability of the data from the two processing strategies, especially taking into consideration their usefulness for climate studies. In the first step of our analysis, the process of screening was conducted to obtain the most homogeneous data sets. Information about various types of oscillations was obtained from the Lomb-Scargle periodograms which were prepared for every station. Then, the value of the seasonal components and linear trends were estimated using the least squares estimation method. Finally, the Mann-Kendall trend test was performed to confirm the presence of a trend in the ZTD time series. On one hand, the spectral analysis showed that the short-term changes of ZTD are almost identical in both campaigns; therefore it can be concluded that the processing strategy does not significantly affect them. Consequently, it can be stated that they can be used as a reliable source of data for climate studies, especially taking into consideration that seasonal oscillations reflect prevailing weather conditions which occur in a given region. Their values reflect the strength of the annual changes in weather, including increasing or decreasing influences of warm water masses, atmospheric pressure centres or other factors affecting the formation of weather. In this sense, due to high spatial resolution, the GPS observations can bring significant benefits to the climate community. On the other hand, long-term changes of ZTD are not so clear. The linear trend in the ZTD time series is particularly important for the purposes of climate studies, because it reflects the increase or decrease of temperature or water vapour content in the atmosphere. The existence of a linear trend in the ZTD time series might, therefore, demonstrate that some physical and real changes take place in the troposphere. Hence, regardless of the selected techniques, its value should be similar or identical. In case of the GPS technique and this study, we are dealing with the use of exactly the same measurement equipment (receivers and antennas). Therefore, the obtained differences could be only caused by the applied processing strategies. In this paper, the linear trend was considered in two lengths of the time series: 16-and 18-years, for both reprocessing campaigns.
Generally, the trends in Repro2 were significantly smaller in comparison to Repro1 for about half of the stations. For the shorter period, the mean trend value in Repro2 was by 0.08 mm yr −1 smaller than in Repro1; however, it is worth noticing that differences in trend values between both campaigns for individual stations were often higher than the linear trend uncertainties. Similarly for the longer period, the difference in mean trend value between Repro1 and Repro2 was 0.09 mm yr −1 and individual differences computed for each station were higher than trend uncertainties. At the same time, the differences between the 18-year and the 16-year ZTD time series were similar in both campaigns. This proves that changing the length of the analysed time series has an influence that is independent of the data sources. The obtained results indicated that the linear changes of ZTD seem to be more sensitive to the applied processing strategy than the seasonal components. The method used for the trend estimation was exactly the same for both campaigns, so the differences in its values resulted only from different approaches to ZTD determination. Investigation of factors which have the biggest influence on the differences is very important in further interpretations. Awareness of the limitations and error sources of ZTD derived from GPS observations, may improve the reliability of the IWV conversion results. Differences in the ZTD values presented in this paper could result from application of various mapping functions which are based on various sources of meteorological data (the Niell mapping functions are based on radiosonde profiles and the Vienna mapping functions 1 are based on the ECMWF weather model). The differences could also be affected by the sequence of determination of the final ZTD value. In Re-pro1, ZTD was fixed to weekly coordinates, while in Repro2 it was estimated together with the other parameters, including coordinates, in daily solutions. Finally, application of atmospheric loadings in Repro2, which have direct impact on the estimated height, could also affect the estimated ZTD. However, it is hard to explicitly separate the influence of each parameter using only the compared data sets, and further analyses are still needed.
Regardless of the differences in the linear trend values we found interesting similarities in the geographical distribution. For both reprocessing campaigns, positive and negative trends occurred in the same geographical area, which indicates that they are not random effects but rather reflect real changes in the troposphere. For the shorter time period, the eastern part of Europe was characterised by positive trends, while in western and south-western Germany, the Netherlands, Belgium and north-eastern France, negative trends were detected. Similar spatial distribution of linear trend characters (between Repro1 and Repro2 campaign) occurred also in case of the longer time period. Moreover, the obtained differences in trend values between the shorter and the longer periods of time indicate that during the analysis of the linear trend for the purpose of climate studies it is necessary to work on exactly the same time span, due to the fact that differences in time series length noticeably affect the value of the trend. This is a particularly important assumption when analysing their spatial distribution.

Data availability
The data sets are available on request.