Statistical significance of trends in Zenith Wet Delay from re-processed GPS solutions

Long series of Zenith Wet Delay (ZWD) obtained as part of a homogeneous re-processing of Global Positioning System solutions constitute a reliable set of data to be assimilated into climate models. The correct stochastic properties, i.e. the noise model of these data, have to be identified to assess the real value of ZWD trend uncertainties since assuming an inappropriate noise model may lead to over- or underestimated error bounds leading to statistically insignificant trends. We present the ZWD time series for 1995–2017 for 120 selected globally distributed stations. The deterministic model in the form of a trend and significant seasonal signals were removed prior to the noise analysis. We examined different stochastic models and compared them to widely assumed white noise (WN). A combination of the autoregressive process of first-order plus WN (AR(1) + WN) was proven to be the preferred stochastic representation of the ZWD time series over the generally assumed white-noise-only approach. We found that for 103 out of 120 considered stations, the AR(1) process contributed to the AR(1) + WN model in more than 50% with noise amplitudes between 9 and 68 mm. As soon as the AR(1) + WN model was employed, 43 trend estimates became statistically insignificant, compared to 5 insignificant trend estimates for a white-noise-only model. We also found that the ZWD trend uncertainty may be underestimated by 5–14 times with median value of 8 using the white-noise-only assumption. Therefore, we recommend that AR(1) + WN model is employed before tropospheric trends are to be determined with the greatest reliability.


Introduction
The position time series derived from Global Positioning System (GPS) observations are widely used to analyze and interpret a wide range of geophysical processes (Teferle et al. 2008;Wöppelmann et al. 2009;Kreemer et al. 2014). Previous studies presented that the noise properties, i.e. the stochastic character, of the GPS position time series are affected by time-correlated noise (Williams et al. 2004; Santamaria-Gomez and Memin 2015; Klos et al. 2016) rather than a pure white noise. It is now widely accepted that the assumption of randomness in the GPS position time series leads to underestimation of the uncertainties of trends, up to an order of magnitude. GPS has proven to be able to infer the conditions in the troposphere to interpret climate and weather (Rohm et al. 2014). A delay in the vertical (zenith) direction known as Zenith Total Delay (ZTD) is estimated from slant total delay of GPS signal. The ZTD has two major components: the Zenith Hydrostatic Delay (ZHD) and the Zenith Wet Delay (ZWD). While ZHD is adequately modeled from surface pressure and temperature, using observations or a model, the ZWD is either directly reported or has to be computed by subtracting ZHD from ZTD. The ZWD can then be converted into the Integrated Water Vapor (IWV) content of the atmosphere (Bevis et al. 1992). Any of these delay estimates have been termed as GPS-derived atmospheric or tropospheric products.
Many applications of the GPS-derived tropospheric products were developed since GPS data contain almost two decades of homogeneously re-processed observations. While initially these data sets were used to analyze meteorological phenomena (Brenot et al. 2006;Labbouz et al. 2013), long-term investigations employing IWV estimates using up to 20 years of GPS solutions became possible (Vey et al. 2009). Nevertheless, the trend and its error estimated considering the stochastic properties of data is still a major source of uncertainty to understand the global climate system. Therefore, we need to correctly characterize the noise in ZWD time series on a global scale as it is crucial to estimate the uncertainties of the trend.
The stochastic properties of climate data are described by an autoregressive noise process (Oladipo 1988;Mann and Lees 1996;Matyasovszky 2012). ZWD is directly linked to climate processes and water vapor variability in the atmosphere, so the noise in ZWD series can arise from measurement and processing errors as well as atmospheric variability. Based on that, we might expect that the measurement and processing errors will introduce a white noise component (Langbein and Johnson 1997), while the atmospheric variability will resemble the autoregressive process (Trenberth 1985). This implies that we can expect that ZWD trend uncertainties will increase compared to the white-noise-only assumption (Combrink et al. 2007;Nilsson and Elgered 2008). This motivates us to assess the stochastic properties of ZWD time series from homogeneously re-processed GPS solutions and deliver reliable estimates of trends and their errors. This analysis will lead to a full understanding of ZWD trends which will be then interpreted as climate changes.
We focus on the stochastic model which is preferred to characterize the ZWD time series. We examine a firstorder autoregressive process and compare it to the pure white noise process used in previous studies. We prove that a combination of an autoregressive process of first order, which has been used by climatologists, with added white noise provides a good description of the power in ZWD data. Finally, we give the statistics of the percentage contribution of the autoregressive process in a combination of autoregressive plus white noise, which clearly indicates the need for autoregressive noise to be applied.
We end with the conclusion of how much one is misled by the improper assumption of noise. This paper is organized as follows. Firstly, we give a description of methodology we used. Also, we focus on ZWD time series we employed. Then, we define the noise models we examined. We describe the results in the 'Results' section. At the end, we summarize, discuss and conclude the main findings we obtained.

Methodology
In the following section, we describe the GPS data processing strategy we have employed to compute the homogeneous daily GPS solutions for this study, including the modeling and estimation of the ZTD and ZWD values. We detail our basic homogenization strategy applied to the ZWD time series and finally, we describe the ZWD noise models we have investigated.

ZWD estimation
As one of International Global Navigation Satellite System Service (IGS) Tide Gauge Benchmark Monitoring (TIGA) Working Groups, the University of Luxembourg completed the newest solution containing up to 750 permanent GPS stations that included all observations from 1995 to 2017. The IGb08 core stations (Rebischung et al. 2012) along with a large amount of GPS stations of TIGA situated near or at tide gauges were included in the processing.
The Bernese GNSS Software 5.2 (BSW52; Dach et al. 2015) was applied to re-process GPS observations in a double difference network strategy with a use of the Centre for Orbit Determination in Europe (CODE) products and the latest International Earth Rotation and Reference Systems Service (IERS) 2010 conventions (Petit and Luzum 2010). We employed the Vienna Mapping Function 1 (VMF1) (Böhm et al. 2006) with coefficients derived every 6 h from the pressure-level data estimated by European Centre for Medium Range Weather Forecasting (ECMWF) (Simmons and Gibson 2000). ZHD was estimated every 1 h using the ECMWF Numerical Weather Prediction (NWP) model. Please see Hunegnaw et al. (2016) for the details of processing. The above processing resulted in more than 20 years of hourly sampled ZWD time series. We employed 120 stations with the length of data from 6 to 22 years, please see online Resource 1 and Fig. 1.

Homogenization of ZWD time series
Although homogenously re-processed, the ZWD time series may be affected by unwanted features such as discontinuities and outliers. Abrupt changes in ZWD may arise from climatic events or be a consequence of GPS data quality issues. For both artifacts, we assume that if the effect is seen in the height position estimate, it is likely also to have propagated into the ZWD estimate.
Homogenization is the technique which combines together several steps of detecting discontinuities, verifying their presence and correcting the climate time series (Ning et al. 2016). It is crucial for ZWD time series to detect all discontinuities as their magnitude is comparable with the magnitude of climate signal and also corrupt trend estimate. In our re-processing, we identified 2500 discontinuities in the 750 position time series based on a manual inspection of all series and used the discontinuity file of the International Terrestrial Reference Frame 2014 (ITRF2014; Altamimi et al. 2016). Two above-mentioned were also supported by epochs estimated with the Sequential t test algorithm (STARS; Rodionov 2004) using a segment length of 100 days and a confidence level of 95%.
We employed a very conservative approach to performing a basic homogenization and used all offsets found for position time series. To reduce the number of epoch candidates, we applied a threshold criteria. For this, we estimated the amplitudes of the ZWD discontinuities with a Least Squares Estimation (LSE) method and only retained those discontinuities with magnitudes higher than three times the estimated offset uncertainty. For simplicity, during this process we assumed a white-noise-only model but investigated any amplitudes suspected to be biased due to numerical artifacts. The third step in our basic homogenization strategy included a manual inspection of the ZWD time series. In this way, the total of 530 offsets reported for the set of 120 stations was reduced to 333 discontinuities that were retained for the correction stage (see Online Resource 2 and Fig. 2).
The basic homogenization strategy also dealt with the detection and removal of outlying ZWD estimates. For this, we applied the Interquartile Range (IQR; Langbein and Bock 2004) approach to remove observations that stand out from others using a threshold of three times the interquartile range. Interquartile range is computed as a difference between 75th percentiles and 25th percentiles or, in other

Functional model of ZWD time series
To derive the trend and amplitudes of seasonal components together with the character of the stochastic part, we employed the reformulated maximum likelihood estimation (MLE) in the Hector software (Bos et al. 2013). We decided on the preferred noise model based on the values of the logarithm of MLE function, which is further referred to as MLE, and Bayesian information criterion (BIC; Schwarz 1978). For the trend and amplitudes of seasonal signals, a time series model can be fitted as: We included k = 6 harmonic terms of frequencies f k = (1/365.25, 2/365.25, 3/365.25, 4/365.25, 1, 1/2) days. The term v is the linear trend, and ZWD R is the initial value at the reference time t i = t R . The coefficients of the harmonic terms are represented by C k and S k . ε ZWD is the stochastic part which is formed when a deterministic part is removed.
First and second harmonics of the annual period were included in the deterministic model for the entire set of stations, as they were proven to be significant in this research using the stacked Power Spectral Density (PSD) computed with the Welch periodogram (Fig. 3). The periods of 121.75 and 91.31 days were assumed only for stations for which their amplitudes were higher than the (1) errors of their estimates. The daily and half-daily curves were assumed for the entire set of analyzed stations, following previous studies by Jin et al. (2008). Other harmonics up to the 12th were also examined, as clearly seen in the spectra, but they were not significant. Their amplitudes ranged between 0.05 and 0.3 mm. Having taken this fact into consideration, we decided not to model harmonics with period shorter than a day. The peaks observed on the individual and stacked PSDs resulted from the effect of the length of the window we employed to compute the Welch periodogram. Having modeled the exact frequencies, we have clearly noticed in the spectra that the power of the peaks did not decrease, which confirmed that they were artificially generated.
The RMS values of ZWD residuals ε ZWD were estimated when mathematical model of series was removed from ZWD observations. It means that trend, offsets, outliers and all seasonalities have been removed prior to RMS estimates. RMS values vary significantly as they depend on the spread of individual series and varied from 1.8 to 50.6 mm, depending on the station.

Stochastic modeling of ZWD time series
We fit different processes into ZWD residuals and compute the BIC and MLE values to decide on a preferred noise Fig. 2 Exemplary ZWD plotted for station TOW2 in Australia. Original ZWD is presented in gray, the deterministic model we employed is plotted in black. We also added the offsets we found in the basic homogenization strategy. These are plotted in red Fig. 3 Stacked spectra of ZWDs for a set of 120 stations. It shows significant oscillations of a frequency of 1, 2, 3 and 4 cycles per year (cpy) as well as 1 and 2 cycles per day (cpd). Although stacked spectra do not show an evident 3-month peak, some individual stations are characterized by this oscillation. The power of the ZWD residuals is plotted in gray. The power averaged with moving window is plotted in red. A set of peaks which may be noticed for high frequencies was found to be artificial according to the length of the window employed to compute the PSD model. We start with the widely assumed pure white noise model (WN). Then, we implement a combination of powerlaw and white noise (PL + WN) processes. In this case, matrix C which is a covariance matrix of observations is re-computed to follow the assumption we made: I and J k are, respectively, the covariance matrices of white and colored noise (Williams 2003). a and b refer to the amplitudes of the white and power-law noise processes.
To follow the climatologists, we also characterize the stochastic part of ZWD data as an autoregressive process (AR). Due to the fact that the stacked spectra of ZWD show a clear flattening in high and low frequencies, we propose to add white noise to the autoregressive process, obtaining a combination of autoregressive and white noise processes (AR + WN).
In general, the autoregressive process of an nth-order AR (n) noise model is defined as: where are the autoregressive (AR) coefficients, while ε means the residual value mentioned in (1). Z i is a Gaussian variable whose standard deviation is known. A ZWD residual time series, ε ZWD (t), follows AR(1) if it is characterized by the relationship (Sowel 1992): where is an autoregressive coefficient of the first order. The AR(1) model refers to the case when any data sample is explained by the previous observation plus a random increment due to additive white noise.
We divided trend uncertainties estimated with AR(1) + WN (σ AR(1)+WN ) and WN (σ WN ) models to deliver the ratios between them: The above quantify shows by how much we were misled when WN was assumed instead of AR(1) + WN.

Results
In the following section, we present results of the deterministic part along with the parameters delivered for the noise analysis. We also compare different noise processes and end with a recommendation on the preferred noise model to be used for any future analysis of Zenith Wet Delay. (

Temporal variations of ZWD
ZWD time series are characterized by different seasonal behaviors due to the area in which the station is located. We examined the stacked and individual PSDs for stations we employed. Figure 3 presents a stacked ZWD PSD for a set of 120 stations along with the averaged power. Clearly noticed is the fact that the hourly sampled ZWD time series contain a significant annual peak and its three overtones. The amplitudes of annual oscillation dominate for the entire set of stations we examined. For 70% of stations, the semi-annual signal is two times smaller than the annual one. Periods of 3 and 4 months, i.e., 4 and 3 cycles per year (cpy), were noticed for low-and mid-latitude stations, while they were hard to be noticed for high-latitude stations. As an example, station NTUS (Singapore, Republic of Singapore) is characterized by a 4-monthly signal which is as powerful as the annual one.
The annual amplitudes varied from 10 to 150 mm (Fig. 4). The largest estimates were obtained for stations in low-and mid-latitudes. The phase shifts were consistent within hemispheres and fall between July and August for the Northern Hemisphere (NH) and between January and February for Southern Hemisphere (SH).
The diurnal cycle of ZWD reflected mainly variations in pressure, temperature, and humidity during the course of the day. For the analyzed ZWD time series, stations situated in low-latitude regions were characterized by the largest diurnal peaks (Fig. 5). These peaks were almost 5-10 times higher than the amplitudes estimated for any other region of the world. The diurnal amplitude for polar stations varied between 0.2 and 1.0 mm. European stations were characterized by a very consistent time of diurnal maxima around 18 h. Stations situated between 30° and 90° for both hemispheres were characterized by similar diurnal amplitude.

Noise analysis of ZWD
The residuals of ZWD time series ε ZWD were analyzed assuming different noise models. The preliminary assessment of a noise model which is preferred for ZWD data was based on the stacked PSDs (Fig. 3). A clear powerlaw behavior was noticed for frequencies between 36.5 days and 8.8 h (10-1000 cpy). For frequencies higher than 1000 cpy and lower than 10 cpy, an evident flattening (or pure WN) was observed. These presumptions were confirmed by the individual spectra. An example is shown in Fig. 6 for two stations situated in different climate zones: MANA (Managua, Nicaragua) and SYOG (Syowa, Antarctica). For residuals computed for any individual station, we fitted four different noise models. We started with a pure WN model, used until now to describe the stochastic part of troposphere data derived from GPS. This model is used as a reference in our research. However, a flat spectrum, which characterizes the WN, is not what ZWD data reveal. The PL + WN model employed for GPS position time series fits medium and high frequencies of the PSD. Failing to capture the lowest frequencies, it leaves some of the power unexplained leading to an artificial increase in the values of uncertainties of the determined parameters. The low-order autoregressive noise models, as supposed at the beginning of this study, are the best for ZWD from all the stations we employed. The use of AR(1) model provides a good characteristic of the frequency spectrum of ZWD residuals. Nevertheless, it is not well fitted to the PSD above 100 cpy. When white noise is super-positioned to AR(1), leading to the combination AR(1) + WN, the best fit of ZWD residuals is provided.
We compared the BIC and MLE values for all noise models we employed. Both values showed that the AR(4) + WN model fitted the ZWD residuals better than WN, PL + WN, AR(n) + WN of orders 1-3 and AR(n) of orders 1-4. However, we have noticed that the differences in BIC and MLE values between certain orders of AR(n) + WN are small and that the statistics themselves are almost equivalent. Based on these results, we present the AR(1) + WN as the preferred noise model for ZWD time series. Going to higher orders of the AR process does not lead to a significant improvement in fitting the spectra, but takes several times as much computation time as AR(1) + WN does, depending on the length of data. Also, we noticed a change in the trend uncertainty of 0.05 mm/ year at maximum. Moreover, we did not find any explanation why higher orders of the AR process should be implemented to model ZWD time series, as climate variability was reported to be the best fitted by low-order autoregressive noise. On the basis of that, we decided to only pursue the AR(1) + WN model further in a comparison to the pure WN model and do not discuss the other models further.
We estimated the fraction of AR(1) noise in the employed combination (Fig. 7). We found that for 103 out of 120 stations, the AR(1) noise model contributes to (1) + WN in more than 50%. This means that it overbalances the white noise which is added to the combination. Also, for 14 stations the percentage of AR(1) contribution is higher than 95%, which means that the ZWD residuals resemble a nearly pure autoregressive behavior. Only for 1 (MCM4, Ron Island, Antarctica) out of 120 stations, the fraction of AR(1) noise is lower than 10%, meaning that WN is dominant. The AR(1) process dominates white noise for stations situated in Europe, New Zealand, Japan and east coasts of Northern and Southern America. WN overbalances AR(1) for few inland Asian stations, Antarctica, the western part of Greenland and a few stations in Northern America. The coefficients ϕ of the autoregressive noise are at the median level of 0.976 ± 0.001 for the entire set of analyzed stations.
The amplitudes of the combination of AR(1) + WN were found to range between 9 and 68 mm with a clear latitude dependence (Fig. 8). Stations with the latitudes higher than 60° are characterized by noise amplitudes smaller than 20 mm. Also, a drop might be observed for stations in the equatorial zone, between latitudes of − 10° and 10°. All other stations are characterized by the amplitudes between 20 and 78 mm. The highest amplitudes (values greater than 60 mm) were noticed for COCO (Cocos Island, Australia), IMBT (Imbituba, Brazil), TOW2 (Cape Ferguson, Australia), MOB1 (Fort Morgan, USA) and J211 (Tsumagoi, Japan). Moreover, the spread of the noise amplitudes is also latitude dependent. A lower spread of amplitudes was noticed for the polar regions for the latitudes beyond 60°, while it was much larger for latitudes between − 60° and 60°. Online Resource 1 lists the trends with associated uncertainties for the WN and AR(1) + WN models. Figures 9 and 10 present an investigation into the significance of trends when the two noise models are employed. The significant trend means here the value of trend higher than its 1−σ error. For a WN assumption, the trends estimated for 115 stations out of 120 were significant. Stations GLSV (Kiev, Ukraine), LAMA (Lamkowko, Poland), LPGS (La Plata, Argentina), RESO (Resolute, Canada) and SAMO (Fagalii, Samoa) were the only exceptions. For AR(1) + WN processes, trends were significant only for a set of 77 stations out of 120, which means that in 43 cases the trend estimates were lower than their 1−σ errors. The ratio between AR(1) + WN and WN trend uncertainties varies from 5 to 14 for stations we employed ( Fig. 11 and Online Resource 1), with an average of 8. Fewer stations' trends became statistically insignificant for the case when we changed the noise model from WN into AR(1) + WN, e.g., for HOFN (Hoefn, Iceland), METS (Metsahovi, Finland), ALIC (Alice Springs, Australia), ZIMM (Zimmerwald, Switzerland). A set of 91 out of 120 used stations has a ratio greater than 8, while 22 stations are characterized by a ratio greater than 10. An uncertainty of trend estimates for ZWD data can be used for climate applications only if the stochastic character, which is not a white noise model, is fully accounted for. An alternative would be to follow the approach presented by Mazzotti et al. (2005) for the position time series. Take the uncertainties of ZWD trends estimated with WN process and multiply them by 8 to estimate the approximate uncertainty that the ZWD trend would have if computed with AR(1) + WN processes. The ratio of uncertainties is strictly related to different factors: the length of residual ZWD time series, the contribution of the AR noise model in the AR(1) + WN combination and, also, the There is no answer to which of those values influences the ratio the most. However, few stations characterized by a high contribution of AR(1) noise to the AR(1) + WN combination show higher ratio of trend uncertainties than any other stations.

Discussion
The time series of ZWD tropospheric data, produced by homogenously re-processed GPS observations, may be assimilated into climate models as they now are into NWP models. The introduction of these atmospheric products showed a measurable improvement in the accuracy of the forecasts (Yan et al. 2009;Mahfouf et al. 2015;Wilgan et al. 2016). Although the assimilation of ZWD data into climate models is a relatively new problem, a perfect example of the process is the ongoing assimilation of the ZWD time series into a model of the Uncertainty Estimation in Regional ReAnalysis (UERRA) Project managed by the UK Met Office.
In this research, we applied the homogenized ZWD time series like the ones used by Jin et al. (2007) to evaluate the parameters of the deterministic and stochastic parts of these ZWD time series. The estimates of amplitudes and phase shifts of seasonal components match the seasonal signals presented in Jin et al. (2007). The daily curve phase shifts also matched the values presented by Ning et al. (2013) with a tolerance of 2 h.
White noise, which has been commonly used to describe the nature of the stochastic part of ZWD, failed to match the power spectral density of ZWD residuals. The PL + WN, which is commonly used to describe the position time series of GPS permanent stations, describes the middle and high frequencies of the ZWD time series spectrum well. However, the PL + WN does not perform well as a noise model to describe the low frequencies of the same spectrum and may cause an artificial correlation between the ZWD time series residuals. The AR(1) model, which is the most often used by climatologists for this application (Percival et al. 2007), conforms very well to the power of ZWD residuals only when complemented by white noise. The AR(1) model matches the medium-frequency spectrum perfectly. White noise in an AR(1) + WN combination results in the model explaining the variability of the spectrum at low and high frequencies. Combrink et al. (2007) demonstrated that the errors of the IWV data trend determined for the HRAO (Hartbeespoort, South Africa) GPS permanent station were significantly reduced when the applied noise model was changed from autoregressive moving average process ARMA(1,1) to white noise. In our research, we used a time span of 19.7 years of data from the HRAO GPS permanent station. The trend error we estimated was 0.01 mm/year for white noise assumption, and 0.13 mm/year for the AR(1) + WN. If a shorter time series were to be used, the error of the determined trend would naturally be much larger. Jin et al. (2007) evaluated the trends for 2-h ZTD data produced from 150 IGS permanent stations and demonstrated that nearly one-half of the determined trends shows that trend errors were underestimated a few times when WN model was assumed prior to analysis. AR(1) + WN model gives much more realistic values of errors of a trend which will be used for climate applications Fig. 11 Ratio of trend uncertainties computed using AR(1) + WN and WN models. For details, we refer to Online Resource 1 exceeded 2 mm/year. We also focused on identifying the significant trend values and showed that nearly all trends determined for the ZWD time series in the Northern Hemisphere were positive, whereas those in the Southern Hemisphere were mostly negative (please see Fig. 9 and online Resource 1). A flaw in these analyses was to apply white noise to attempt to describe the nature of the stochastic part of the data.
We evaluated the trend errors for a set of 120 permanent stations by applying the AR(1) + WN model. A comparison between the errors we determined and those from Jin et al. (2007) reveals that the latter were underestimated by up to 14-fold. By applying a preferred noise model, i.e. AR(1) + WN, instead of the white noise model, we showed that some trends became insignificant because their values were lower than the corresponding error. Because of this, these trends should not be considered in climate analyses or interpretations. With the nature of the stochastic part of the data assumed to be WN, only 5 out of 120 trends analyzed were statistically insignificant.
The trends determined for the permanent stations NANO, MDO1, ALGO, HOFN, LAMA, KIR0, and WSRT are considered to be insignificant with AR(1) + WN applied. These were also determined using white noise and elaborated by Jin et al. (2007). Some of the permanent stations located in and near Antarctica, e.g., CAS1, and demonstrated earlier by Thomas et al. (2011), had their trends rendered insignificant when the AR(1) + WN model was applied to describe the stochastic part of the data.

Conclusions
In this research, we provided a deep analysis of the character of the ZWD residuals from a set of 120 GPS permanent stations located around the globe, which is required prior to their assimilation into climate models. The ZWD time series applied herein were produced by re-processing the observations gathered from a global network over the period 1995-2017. The application of uniform models to the entire GPS observation span left the resultant ZWD time series without discontinuities that are associated with potential changes or improvements in models over the years (Steigenberger et al. 2007). We performed a basic homogenization of the ZWD time series to detect any other discontinuities arising from changes in the antennae or software through years. This task is critical, as failure to properly remove discontinuities or to complete data homogenization of the ZWD time series results in the incorrect evaluation of the ZWD time series trends and their errors.
To date, the stochastic part of the ZWD time series has been modeled as a pure white noise. It represents the residuals, observed in successive epochs, are uncorrelated with one another in a time domain. White noise has a negligible effect on the errors in the determined parameters of the ZWD time series, including the trend. This is why we propose the model which combines an autoregressive process with white noise. We show that this combination is a preferred choice to describe the nature of the stochastic part of the ZWD time series. The analytical result we produced matches the results published by climatologists whose work is focused on the analysis of different climate data, including air temperature, pressure and/or humidity, by describing them explicitly as an autoregressive process. Henderson-Seelers and McGuffie (2012) noted that a loworder autoregressive process is a preferred one with which to describe atmospheric data. Similar to climate data, the residuals of the ZWD time series also revealed a correlation between the successive observations. The correlation may be modeled very well by combining white noise with an autoregressive process. The model thus formed may be applied with success to all climate zones. For most of the GPS permanent stations we employed, the contribution of the autoregressive process in the modeling combination exceeded 50%, which may have been caused by variability in the atmosphere. The remainder of the residuals were explained by white noise, introduced into the ZWD time series during the gathering of GPS observations or during the processing of the latter.
Our research allows us to conclude that a combination of autoregressive process plus white noise is a model well capable of describing the nature of the stochastic part of ZWD tropospheric time series and should always be assumed in all determinations of these trends. If the assumption about the nature of the ZWD time series is incorrect, the results produced should not be analyzed for any climatic application, since the resulting error of the ZWD time series may be underestimated fivefold to 14-fold when compared with the ubiquitously used white noise assumption. The method we presented can be employed to evaluate tropospheric trends from any other techniques than GPS that belong to GNSS or from any other space technique. mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.