The uncertainty of the atmospheric integrated water vapour estimated from GNSS observations

Within the Global Climate Observing System (GCOS) Reference Upper-Air Network (GRUAN) there is a need for an assessment of the uncertainty in the integrated water vapour (IWV) in the atmosphere estimated from ground-based global navigation satellite system (GNSS) observations. All relevant error sources in GNSS-derived IWV are therefore essential to be investigated. We present two approaches, a statistical and a theoretical analysis, for the assessment of the uncertainty of the IWV. The method is valuable for all applications of GNSS IWV data in atmospheric research and weather forecast. It will be implemented to the GNSS IWV data stream for GRUAN in order to assign a specific uncertainty to each data point. In addition, specific recommendations are made to GRUAN on hardware, software, and data processing practices to minimise the IWV uncertainty. By combining the uncertainties associated with the input variables in the estimations of the IWV, we calculated the IWV uncertainties for several GRUAN sites with different weather conditions. The results show a similar relative importance of all uncertainty contributions where the uncertainties in the zenith total delay (ZTD) dominate the error budget of the IWV, contributing over 75 % of the total IWV uncertainty. The impact of the uncertainty associated with the conversion factor between the IWV and the zenith wet delay (ZWD) is proportional to the amount of water vapour and increases slightly for moist weather conditions. The GRUAN GNSS IWV uncertainty data will provide a quantified confidence to be used for the validation of other measurement techniques.


Introduction
In the hydrological cycle, water vapour is an important variable for transferring heat energy from the Earth's surface to its atmosphere and in moving heat around the Earth.Meanwhile, water vapour is a very important greenhouse gas due to its ability to absorb long-wave thermal radiation emitted from the Earth's surface.Hence, the atmospheric water vapour is very important for the Earth's climate system, and its variability is a key to understanding the hydrological cycle.A variety of systems exist for measuring the atmospheric integrated water vapour (IWV), e.g.radiosondes (Ross and Elliott, 1996), microwave radiometry (Elgered, 1993), very long baseline interferometry (VLBI) (Heinkelmann et al., 2007), and global navigation satellite systems (GNSSs) (Bevis et al., 1992).Standing out from others, GNSS observations from the ground can be made under in principle all weather conditions with a high temporal resolution (a few minutes), and its spatial resolution continuously improves when more continuously operating stations are installed.Therefore, using ground-based GNSS measurements to de-Published by Copernicus Publications on behalf of the European Geosciences Union.
T. Ning et al.: The uncertainty of the GNSS-derived IWV rive IWV has drawn a lot of scientific attention and has been investigated (e.g.Wang et al., 2007;Vey et al., 2010;Ning et al., 2013).In order to interpret GNSS measurements correctly and draw valid conclusions on the quality of the resulting IWV estimates, the uncertainty of the GNSS-derived IWV must be carefully evaluated and quantified.
The procedure for the estimation of the atmospheric IWV from GNSS measurements begins with the refraction of the GNSS signals (e.g.GPS signals transmitted at frequencies L1 = 1.575GHz and L2 = 1.228GHz) when passing through the Earth's neutral atmosphere.The refraction introduces an additional delay to the primary observable, the signal propagation time.The propagation delay can be estimated in the GNSS data processing, together with other parameters such as coordinates, as a zenith total delay (ZTD) which is normally expressed in units of length using the speed of the light in vacuum.ZTD is separated into two parts: the zenith hydrostatic delay (ZHD) and the zenith wet delay (ZWD).The ZHD can be determined from surface pressure measurements (Davis et al., 1985), and the ZWD depends on the amount of water vapour in the column of air, i.e. the IWV, through which the signal passes and can be estimated from the data themselves.
All GNSS measurements are subject to error sources that influence the uncertainty of the estimated ZTD which are converted from the slant delays using mapping functions (MFs; Boehm et al., 2006).Since the slant delay is referring to the delay in the slant range distance between the GNSS receiver and a satellite, the satellite orbit errors will propagate to the slant delay and the resultant ZTD.The accuracy of the ZHD depends on the accuracy of the ground pressure.Additionally, the conversion from the ZWD to the IWV will add uncertainties, e.g.due to imperfect determination of the mean temperature of the atmosphere.All those uncertainties shall be included when we calculate the final uncertainty of the GNSS-derived IWV.Errors in GNSS measurements can be random or systematic, or more commonly a mixture of both, depending on the timescale studied.Since the expected (mean) value of random errors is zero, the impact of such errors is reduced as the number of measurements increases.Systematic errors cannot be averaged out as the time series becomes longer.They can however change at a specific time epoch.For example, a change of the GNSS antenna or its environment may introduce such an offset.Systematic errors may also change slowly with time; e.g. the impact of signal multipath at a GNSS site may vary due to growing vegetation (Pierdicca et al., 2014).
Depending on different applications, the requirement on the accuracy of the estimated IWV varies.For the forecasting application the demand is mainly the timing of moving air masses, while the accuracy of individual IWV estimate and the IWV biases is of less importance.Therefore, less accurate real-time orbit estimation is accepted in forecasting applications.For climate research, it is crucial to have a high accuracy with the smallest biases possible in order to obtain the absolute value over a long timescale.In this case, the final GNSS orbit estimation with the highest accuracy is necessary.The focus of this study is to discuss and assess the accuracy of the IWV derived from ground-based GNSS measurements obtained from post-processed data and mainly used for climate research.
Many studies have investigated the uncertainty of the GNSS-derived IWV either by comparing the GNSS IWV with data obtained from other independent techniques (e.g.Pacione et al., 2002;Rózsa, 2014) or by combining the uncertainties of the input variables used in order to obtain the uncertainty of the GNSS IWV (e.g.Wang et al., 2007;Ning et al. (2013);and Van Malderen et al., 2014).In order to conduct a complete and comprehensive investigation on the uncertainty of the GNSS IWV, the two methods are both discussed in this study in order to develop a method which is applicable to each individual data point.The uncertainty method will be implemented to the GNSS IWV data stream for the GCOS (Global Climate Observing System) Reference Upper-Air Network (GRUAN).The objectives of GRUAN are to provide long-term, high-quality upper-air climate records, to constrain/calibrate data from more spatially comprehensive global observing systems (including satellites) and to measure a large suite of co-related climate variables (Bodeker et al., 2015).One of the most important objectives of GRUAN is to collect reference observations.Each GRUAN reference observation is required to provide a comprehensive uncertainty analysis (Immler et al., 2010).Ground-based GNSS IWV was identified as a priority 1 measurement for GRUAN.This study is intended to develop methods to calculate GNSS IWV uncertainties for each data point and to work with the GRUAN GNSS data processing centre, the German Research Centre for Geosciences (GFZ), to incorporate them into the GRUAN GNSS data stream.
The paper is structured as follows.Sections 2 and 3 discuss the two uncertainty analyses: a statistical and a theoretical one.A sequence of subsections is given in Sect. 3 in order to describe each of the errors which contributes to the final total uncertainty of the GNSS-derived IWV.In the last subsection, the individual error sources are then summarised into an overall error budget in order to obtain the final uncertainty of the GNSS-derived IWV.The conclusions are given in Sect. 4.

Statistical analysis
A statistical analysis (Moses, 1986) can be used to evaluate the uncertainty of the GNSS-derived IWV if independent estimates are available from at least three co-located techniques, measuring the same true variability of the IWV.The individual IWV estimates from techniques A and B are expressed as Table 1.ZWD and corresponding IWV uncertainties given by a statistical analysis using observations from Onsala on the Swedish west coast.

Standard deviation a
[mm] VLBI absolute offset +2.The uncertainty of the IWV given by dividing σ ZWD by 6.5 (a mean value of the conversion factor Q, given by radiosonde measurements, for the Swedish west coast).The uncertainty of Q was neglected.f The uncertainty of the IWV calculated in a similar way as σ V1 but taking the uncertainty of 0.1 kg m −2 for Q into account (σ V2 = σ 2 V1 + 0.1 2 ).The uncertainty of the ground pressure is insignificant since the ground measurements, accurate to ±0.1 hPa, were used.g Assumed values in order to relate the comparisons to an absolute scale. (1) where V i is the true value of the IWV for the epoch i; M A and M B are the biases (the systematic error) for each technique, respectively; and ε A(i) and ε B(i) are the random errors.The observed standard deviation (SD) of the IWV difference between techniques A and B based on N simultaneously sampled data points can be expressed as where the overline denotes the mean value.From Eqs. ( 1) and (2) we have where V is the mean value of the true IWV.If we combin Eqs.
(1), (2), (4), and (5) with Eq. (3), we get For a long time period of measurements giving a zero mean of random errors (ε A = ε B = 0), and assuming that ε A and ε B are uncorrelated, Eq. ( 6) can be expressed as After including the third technique C, we will have two more equations for SD: where the ε terms can be solved using Eqs.( 7) to (9).Since the ε terms only address the random errors, in order to determine the total uncertainty of the IWV, caused by both random and systematic errors, the biases (M A and M B ) for two techniques also need to be derived.The mean IWV difference between two techniques is expressed as where the output of Eq. ( 10) is highly dependent on the length of the timescale of measurements studied.For a short

T. Ning et al.:
The uncertainty of the GNSS-derived IWV time period where the length of measurements is not long enough to average out random errors, Eq. ( 10) gives a mixture of random and systematic errors.For a long time period giving a zero mean of random errors, the equation gives M A − M B .A typical formal error of the GPS ZTD is around 3 mm, which approximately corresponds to a formal error of 0.5 kg m −2 in IWV.Assuming that we have one data point for each day, we need at least 2 years of data in order to reduce the formal error below 0.001 kg m −2 .From Eq. (10), for a long time period of measurements, the IWV bias (M) for one technique can only be obtained when it is compared to a technique with a known IWV bias or no bias at all.After obtaining the values for both ε and M, the total uncertainty of the IWV can be determined: One example of a statistical analysis can be found in Ning et al. (2012).The authors presented the results from comparisons of 10-year-long time series of ZWD, estimated using GNSS, geodetic VLBI, and a water vapour radiometer (WVR).These three techniques are co-located at the Onsala Space Observatory separated by less than 100 m.Since the VLBI antenna has a very high directivity, the VLBI-derived ZWD is not affected by signal multipath.Meanwhile, no instrument changes have occurred on the Onsala VLBI antenna for the last 30 years.Therefore, the VLBIderived ZWD for the 10-year-long time period can be considered as having only a small constant bias.Using Eqs. ( 7) to (11), we calculated total uncertainties in the ZWD and in the corresponding IWV, for each technique after assuming a bias of +2, 0, and −2 mm in the VLBI-derived ZWD.The results are summarised in Table 1, where, depending on the assumed VLBI bias, the uncertainties for the GNSS-derived IWV are around 0.5, 0.7, and 1 kg m −2 , respectively.
Since all known biases must be removed before derivation of the uncertainty, an assumption is made about no bias in another measurement technique.For most of the GRAUN sites, measurements simultaneously acquired from at least three co-located independent techniques are difficult.In this case, the next approach to calculate GNSS-derived IWV uncertainty is more useful.

Theoretical analysis
A theoretical analysis is desired where the total uncertainty of the GNSS-derived IWV (σ V ) is calculated from the uncertainties associated with the input variables according to the rule of uncertainty propagation for uncorrelated input quantities.The equation used for a theoretical analysis is (Immler et al., 2010) where f (v 1 , . . ., v N ) is the functional relationship between the GNSS-derived IWV and the input variables; σ i is the one sigma uncertainty of the corresponding variable.
The GNSS-derived IWV, denoted by V , is converted from the ZWD, denoted by z w , via the conversion factor Q: where the ZWD is given by the subtraction of the ZHD, denoted by z h , from the ZTD, denoted by z t : Combination of Eqs. ( 13) and ( 14) gives the functional relationship Insertion of Eq. ( 15) into Eq.( 12) gives the expression for the propagation of the uncertainties from different sources to the total IWV uncertainty: In order to calculate σ V , we first need to study and determine values for Q, σ ZTD , σ ZHD , and σ Q .

Error budget of the GNSS-derived ZTD
The fundamental observable of the GNSS technique is the propagation time of a signal, transmitted from the GNSS satellites, passing through the Earth's atmosphere to the receivers (ground-based).When we consider the error budget of the GNSS-derived ZTD, errors from several parts need to be taken into account, which will be discussed in this section.

GNSS satellite orbits
Errors in the estimates of the satellite coordinates will propagate directly to the estimates of the GNSS parameters.If we use the precise point positioning (PPP) strategy to process the data obtained from a permanent GNSS site where the site coordinates are usually kept fixed (one estimate per day), eliminate the ionosphere delay to the first order, and use the final clock product from the International GNSS Service (IGS), the orbit error is compensated by the ZTD, receiver clock, and ambiguity parameters (Douša, 2010): where z i A is the zenith angle seen from the ground receiver A to the satellite i; R A and R i are the distances from the top of the receiver to the geocentre and the geometrical distance from the satellite to the geocentre, respectively; σ X i Rad and σ X i Tan are the radial and the tangential orbit errors, where the tangential orbit error is expressed by two vectors: along-track and cross-track; σ ZTD is the ZTD error; σ T A is the receiver clock error; σ N i is the ambiguity error; and λ f is the wavelength of the GNSS signal.Based on Eq. ( 17), Fig. 1 depicts the impact of different components of the orbit errors on the estimated ZTD as the function of zenith angle, and the impact factor indicates how the ZTD error correlates with the orbit error.The factor was calculated using Eqs.( 17) and ( 18) and took the unit (1) for the orbit error.It is evident that the radial orbit error dominates the resulting ZTD error and that the maximum impact occurs when the satellite is in the zenith direction.The maximum impact factor due to the tangential orbit errors is 0.13 and is observed at a zenith angle of 45 • .
Figure 1 was produced by ignoring the parts for the receiver clock error and the ambiguity error in Eq. ( 17) and by assuming that the estimated ZTD compensates for all the radial and tangential orbit errors.In order to determine how many of the orbit errors will be mapped into the estimated ZTD when the two other parameters (receiver clock and ambiguity) are estimated, we carried out simulations for three GRUAN sites.These are LDB0 (Lindenberg, 14.1 • E, 52.2 • N), LDRZ (Lauder, 169.7 • E, −45.0 • N), and NYA2 (Ny-Alesund, 11.9 • E, 78.9 • N).We first processed the GPS data for the year 2014 using a PPP strategy.For each site, we calculated satellite-receiver range based on the estimated site coordinate and the satellite coordinates provided by IGS for each epoch where the corresponding estimated ZTD and receiver clock error from the PPP processing are added back to the calculated range together with known ambiguities.Finally, for each satellite-receiver range we added the range error projected from the corresponding orbit error which was calculated by the left part of Eq. ( 17).The ZTD, receiver clock, and the ambiguity parameters were estimated again using PPP and compared to the old PPP solution in order to determine how the orbit errors map to the different estimated parameters.
Table 2, taken from http://igscb.jpl.nasa.gov/mail/igsmail/2010/msg00001.html,presents the orbit accuracy for each component using the IGS reprocessed orbit product, which should be comparable to the operational final orbits.The mean and standard deviation were calculated based on satellite position repeatability for each pair of consecutive days.For the simulation, we used constants 15 and 38 mm for the radial and the tangential (square root of along and cross) orbit errors, respectively.
Figure 2 depicts results for the three GRUAN sites where the ZTD errors from the simulation for one day (1 June 2014) are plotted along with two other theoretical ZTD errors.One is calculated by Eq. ( 17) but only taking the radial orbit errors   ambiguity).This is due to the fact that major part of the orbit errors in the satellite-receiver range is constant and is therefore to a large extent absorbed by the estimated ambiguities and the receiver clock parameters (Douša, 2010).The mean simulated ZTD errors vary from 1.5 to 3 mm, which show a higher correlation (∼ 0.60) to the theoretical ZTD errors caused by the radial orbit errors than the one caused by the tangential orbit errors.The larger uncertainties, seen for the site LDB0 and after hour 15, are caused both by the larger radial orbit errors and by a worse geometry of the satellites (smaller number of satellites visible).Similar results are seen (not shown) if we perform the simulation for 1 year of data for the three sites.The mean ZTD simulated errors are around 3 mm for the sites LDB0 and LDRZ, while the one for the site NYA2 is 1.5 mm.The corresponding standard deviations are around 1 and 0.4 mm, respectively.

Ionospheric delay
The Earth's ionosphere contains electrons in sufficient quantity to significantly delay the propagation of GNSS signals.The ionospheric delay is dependent on the total amount of free electrons along the propagation path, named total electron content (TEC), and on the carrier frequency of GNSS signals.Normally, in order to remove the ionospheric impact, an ionosphere-free linear combination is used: where 1 and 2 are carrier phase measurements from two different frequencies (e.g.L1 and L2 for GPS).Since this virtual measurement can only eliminate the ionospheric delay to the first order (∼ 99.9 % of the total delay), smaller contributions from higher-order terms remain, which how-ever have more dramatic impact during strong solar activities, e.g.ionospheric storms (Pireaux et al., 2010).Fritsche et al. (2005) carried out an investigation on the impact of including corrections of the ionospheric delay for higher-order terms on GNSS estimates using a global network, and with a focus on the solar maximum between 2001 and 2003.The authors found a linear dependency between the difference in the estimated vertical position and the peak electron densities (varying approximately from 2 to 12 mm when the daily means of the TEC unit increase from 25 to 175) when no corrections were applied for second-and third-order terms.The ZTD difference is approximately a factor of 3 smaller than the vertical position difference (Hill et al., 2009).The corresponding ZTD difference is then from 0.6 to 4 mm.This is consistent with the result reported by Fritsche et al. (2005), where a ∼ 2 mm difference in the ZTD was observed when the TEC unit (a TEC unit equals 10 16 electons m −2 ) reaches 150.However, over long time periods the impact on average values is less.A study was performed by Petrie et al. (2010) for the time period from 1995 to 2008.Using approximately 60 globally distributed GNSS sites, they found that the mean difference in the vertical component, after applying second-and third-order terms, was in the range from −0.3 to 0.5 mm (corresponding to the ZTD difference from −0.1 to 0.2 mm) for sites with data spanning at least 10 years.The result indicates that the impact on the ZTD from the higher-order terms is insignificant over a long time series.

Signal multipath
Differences in the ZTD caused by signal multipath depend strongly on the elevation angle of the observation and are different from site to site due to a changing electromagnetic environment.In order to demonstrate the impact of multipath effects, we carried out analyses using different elevation cutoff angles in the PPP data processing for three GRUAN sites (LDB0, LDRZ, and NYA2) together with two other sites: LDB2 (14.1 • E, 52.2 • N) and POTS (13.1 • E, 52.4 • N).One month of GPS data acquired from June 2014 were processed using four different elevation cutoff angles: 5, 10, 15, and 20 • .We calculated the monthly mean of the estimated ZTD for each solution relative to the value from the solution using a 5 • cutoff angle.The results are shown in Fig. 3, where the most significant elevation-dependent ZTD difference is seen for the site LDB0, while the least is given by the site LDRZ.The results indicate that the impact of multipath effects is highly site-dependent.This is further examined by the two co-located sites: LDB0 and LDB2, which are both located at Richard Aßmann Observatory of the German Weather Service, Lindenberg, Germany, with a baseline less than 100 m.The photographs of the two sites are given in Fig. 4. The two sites use two different antennas.The one used for LDB2 is LEICA25, which is a choke ring antenna, while the one used for LDB0 is a JAV_GRANT_G3T (www.ngs.noaa.gov/ANTCAL/main.jsp).The choke ring antennas are designed to reduce multipath signals.The result for LDB2 shows fewer multipath effects than the one for LDB0 (see Fig. 4).However, the vertical metal rod beside the LDB2 may cause multipath, while the flat surface below the antenna may also cause systematic effects.Furthermore, the LDB2 antenna is covered by a plastic radome, which may also introduce elevation-dependent impacts (see discussion in next section).
A similar investigation on the multipath effects was carried out by Ning et al. (2011) using GNSS sites at the Onsala Space Observatory.Using nearly 80 days of GNSS data, they found an IWV difference of 1.7 mm (corresponding to ∼ 11 mm in ZTD) between a 5 • solution and a 20 • solution.Such a difference can be significantly reduced by using microwave-absorbing material, e.g.ECCOSORB ® , attached below the antenna plane.
Signal multipath is highly dependent on the local environment and can vary in time, e.g.due to changes in soil moisture (Larson et al., 2010) and growing vegetation (Pierdicca et al., 2014).Therefore, it is difficult to have a general model to quantify the multipath effect on the estimated ZTD.For all GRUAN GNSS sites, a microwave-absorbing material, such as ECCOSORB, is highly recommended to be installed to the antenna in order to minimise the impact of the multipath effects.In addition, a careful and continuous documentation of each site is necessary not only for the change of antennas, radomes, and receivers but also for changes in the environment close to the antenna, e.g. a growing vegetation and/or the cutting of trees.

Antenna-related errors
In order to obtain the highest accuracy in the ZTD estimates, antenna-related errors, i.e. phase centre variations (PCVs) and radome effects, need to be removed.Therefore, an absolute calibration of the PCV for all the GNSS satellitetransmitting antennas and the ground antenna (Schmid et al.,  2007) is necessary to be implemented in the GNSS data processing.The difference in the ZTD estimated with and without applying the PCV correction can vary from 2 to 10 mm (Byun and Bar-Sever, 2009;Thomas et al., 2011).Nowadays PCV corrections are normally always applied in order to eliminate this error source.
To avoid the accumulation of snow and for general protection, many GNSS antennas are equipped with radomes.Different shapes of radomes yield different impacts on the phase of the GNSS signal.Emardson et al. (2000) found that the IWV offset introduced by a conically shaped radome can be up to 1 mm (about 6.5 mm in ZTD).A smaller impact (less than 2 mm in ZTD) was observed by Ning et al. (2011) from the use of a hemispheric radome, which is recommended if a GRUAN GNSS site is going to install a radome.Ning et al. (2011) also found that the impact of the radome depends on the different geometries of the electromagnetic environment of the antenna, as well as on the elevation cutoff angle for the observations used in the analysis.Therefore, applying a radome calibration in the data processing is necessary in order to reduce the radome impacts.This is done by using calibration tables provided and updated by IGS which include a set of calibrations for a particular antenna and for a particular radome.

Mapping functions
In GNSS data processing, the slant path delay is converted to the equivalent ZTD (sum of the ZHD and the ZWD) using MFs: where is the elevation angle seen from the ground receiver to the satellite; m h and m w are the hydrostatic and the wet mapping functions, respectively.21) and ( 22) were taken from Table 1 in Stoew et al. (2007).
tions (NMFs) (Niell, 1996) and the improved mapping functions (IMFs) (Niell, 2000) using radiosonde data.The resulting mean (M) and standard deviation (σ ) of the slant delay errors were modelled by Values calculated from Eqs. ( 21) and ( 22) for four elevation angles (5; 7; 10; 15 • ) are summarised in Table 3.It is clear that the accuracy of MFs is highly elevation-dependent and that the MF-induced errors on slant delays are insignificant when the elevation angle gets larger.
Currently the most popular and accurate MF is the Vienna Mapping Function 1 (VMF1) since it can capture the shortterm variability of the atmosphere by utilising data from a numerical weather model (Boehm et al., 2006).Zus et al. (2015) carried out an investigation on the accuracy of the VMF1 by comparing a MF called GFZ-VMF1, which is based on the VMF1 concept, with a direct mapping approach.They found that, if a GPS site is located 2 km above the MF model orography and if a low elevation is used (3 • ), the error on the estimated vertical component can be up to 1 cm (corresponding to approximately 3 mm in ZTD).They also pointed out that it is difficult to distinguish the MF-caused error from a variety of other errors presented at the low elevation angles, e.g.poor or missing antenna PCV models and multipath.Given this fact, for the GRUAN data products, instead of quantifying the ZTD uncertainty from the applied MF, we may choose a slightly higher elevation cutoff angle, i.e. > 10 • , for the GPS data processing in order to significantly reduce the MFinduced error.

Summary of the ZTD uncertainty
All GRUAN GNSS data will be processed by GFZ using its inhouse GNSS software package, Earth Parameter and Orbit determination System (EPOS) (Deng, 2012), where only the PPP approach will be implemented in the data processing.The formal error, provided for each time epoch by the estimation process of PPP, is only dependent on the amount of carrier phase measurements and the constellation of the satellites for a given site (Byun and Bar-Sever, 2009).In order to take systematic errors in the GPS orbit into account, the method discussed in Sect.3.1.1will be applied and the calculated ZTD error for each time epoch will be added to the corresponding formal error.In addition, in order to reduce the ZTD uncertainty due to other factors, i.e. the ones discussed in Sects.3.1.2-3.1.5,the following conditions need to be fulfilled in the GPS data processing: -Corrections for the second order of the ionospheric delay are applied.
-Final orbit/clock products from IGS or equivalent are used.
-Absolute satellite and ground antenna PCV models and radome calibrations are implemented.A hemispheric radome and an ECCOSORB plate are recommended.
-Signal multipath effects need to be minimised either by implementing microwave-absorbing material to the antenna or by locating the GNSS antenna in an favourable place.
-Use an elevation cutoff angle of 10 • or higher.A higher elevation cutoff angle will degrade the geometry and increase the formal error of the ZTD estimate.However, this may still be desired for applications where longterm trends are estimated and systematic errors such as signal multipath rather than formal errors are the limiting factor (Ning and Elgered, 2012).
The procedure to determine the total ZTD uncertainty for each time epoch is described as follows.For each GRUAN site, the daily GPS data will be processed first using a PPP strategy in order to obtain the ZTD estimates and corresponding formal errors together with receiver clock errors.Thereafter, the estimated ZTD and receiver clock errors will be used in a simulation in order to estimate the ZTD errors due to the orbit errors.Then the root sum square (RSS) of the simulated ZTD error and the corresponding formal error gives total ZTD uncertainty for each time epoch.Figure 5 depicts the estimated ZTD for one day (1 June 2014) for three GRUAN sites together with the corresponding calculated total ZTD uncertainty for each time epoch (5 min).The results show that the daily mean ZTD uncertainty is 3.8 mm for the site LDB0, while the values for the sites LDRZ and NYA2 are 3.7 and 3.0 mm, respectively.The result also indicates that the ZTD uncertainty is not correlated with the magnitude of the ZTD.These values are consistent with the claimed 1σ uncertainty (4 mm) by IGS's final ZTD product.

Uncertainty of the ZHD
The ZHD for a given GNSS site can be calculated using the ground pressure z h = (2.2767± 0.0015)  where the constant 2.2767 is calculated from the refractivity of dry air, the universal gas constant, and an effective value of the acceleration due to gravity at the site (Bevis et al., 1994).
The resulting uncertainty includes the different contributions which are assumed to be uncorrelated.In Eq. ( 23), ZHD is in units of millimetres (mm), P 0 in hectopascals (hPa), and where λ and H are the site latitude in degrees and the height above the geoid in metres, respectively.The expression f (λ, H ) is used to model the variation of the acceleration due to gravity in the ZHD.An example, for different latitudes and heights, is shown in Fig. 6.The uncertainties in the latitude and the height are negligible.The uncertainty of the ZHD is mainly determined by the uncertainty in the ground pressure (σ P 0 ) and the uncertainty of the constant (σ c ).Using Eq. ( 12), we have the uncertainty of the ZHD in millimetres: All GRUAN sites should be equipped with surface barometers which provide accuracies much better than 0.5 hPa.One example is shown in Fig. 7, which depicts the pressure differences between three barometers at the Onsala Space Observatory, Sweden.The barometer from the Swedish Meteorological and Hydrological Institute (SMHI) was used as a reference.The absolute accuracy of the SMHI barometer is T. Ning et al.: The uncertainty of the GNSS-derived IWV better than 0.2 hPa.It has been calibrated approximately every second year.Additionally the SMHI barometer is traceable to the SI unit within 0.1 hPa.The small annual variations observed from the Setra barometer are likely due to a sensor sensitivity to the temperature.The uncertainty of surface pressure measured by GRUAN barometers will be estimated using uncertainty estimation protocols presented in Immler et al. (2010).It will be used to estimate the ZHD uncertainty using Eq. ( 25).
For the GNSS sites which are not equipped with barometers, other methods have been investigated.Wang et al. (2007) found an average root mean square (rms) difference of 1.7 hPa given by comparisons between the ground pressure interpolated, both spatially and temporally, from nearby surface synoptic observations and local ground measurements at 48 globally distributed GNSS sites covering an 8-year time period.Similar comparisons were carried out by Heise et al. (2009), but using the interpolated ground pressure from the ECMWF meteorological analysis at more than 60 global GNSS sites.Using 1 year of data, they obtained a better agreement with an overall mean bias and a standard deviation of 0.0 and 0.9 hPa, respectively.We did a similar test, but only for the GNSS site at the Onsala Space Observatory, using over 14 years of data.The result shows a mean bias and a standard deviation of 0.1 and 0.6 hPa, respectively.It should be pointed out here that, if we use the ground pressure obtained from numerical models interpolated to the height of the GNSS site, the uncertainty of the GNSS height will have an impact on the derived ZHD.Snajdrova et al. (2005) found that 10 m of height difference approximately causes a difference of 3 mm in the ZHD.

Uncertainty of the conversion factor Q
The conversion factor Q is defined by where ρ w is the density of liquid water, in units of kg m −3 ; R w is the specific gas constant for water vapour, in units of J kg −1 K −1 ); the values of k 3 and k 2 can be estimated from laboratory experiments, in units of K hPa −1 and K 2 hPa −1 , respectively; and T m is a mean temperature, in units of K.In Eq. ( 26), the uncertainties of the density of liquid water (ρ w ) and the specific gas constant for water vapour (R w ) are 0.002 kg m −3 (Wolf, 2008) and 0.008 J kg −1 K −1 (http://physics.nist.gov/cuu/Constants/index.html),respectively.Since the impact of the uncertainties from these two parameters is insignificant (less than 0.1 % of the total Q uncertainty), the uncertainty of Q is determined by the uncertainties in T m , k 3 , and k 2 .The combination of Equation ( 26) and Eq. ( 12) gives In order to evaluate the impact of the uncertainty in each variable on the total uncertainty of Q, we need to specify values and uncertainties for the constants k 3 and k 2 , which are taken from Bevis et al. (1994).Figure 8 depicts the uncertainty of Q as a function of the uncertainty in T m assuming that the value of T m is 279 K. Three groups of values were used for the uncertainties of the constants k 3 and k 2 : the nominal value (given by Table 1 in Bevis et al., 1994), the double of the nominal value and half of the nominal value.It is evident that, when the uncertainty of T m is sufficiently large, it tends to dominate the uncertainty of Q.Furthermore, the uncertainty of k 3 has a larger impact compared to the uncertainty of k 2 .
The parameter T m can be calculated from the vertical profile of water vapour pressure (p w ) and the physical temperature (T ): Globally, T m can be derived from numerical weather prediction (NWP) models.An rms difference of 1.3 K in T m was claimed by Wang et al. (2005) based on global comparisons between the NCEP/NCAR reanalysis and the radiosonde measurements using 6 years of data (NCEP: National Centers for Environmental Prediction of the US weather service; NCAR: National Center for Atmospheric Research).The authors also found a better agreement (1.1 K) in T m obtained from the ECMWF 40-year reanalysis, ERA-40 (Uppala et al., 2005).Due to the fact that the ECMWF data provides a temporal resolution of 6 h and a horizontal resolution of about 50 km and there is normally a difference between the model height and the GPS height, a temporal, horizontal, and vertical interpolation of the ECMWF data to the time and position of the GPS site is necessary.Details for the interpolation of the ECMWF data can be found in Heise et al. (2009), which are summarised as follows.For the horizontal interpolation the ECMWF interpolation library (EMOSLIB, http://www.ecmwf.int)was used.The 6-hourly ECMWF data were linearly interpolated to the temporal resolution of the GPS ZTD (1 h).The strategy used for the vertical interpolation depends on whether the GPS height is above or below the lowest ECMWF level.For the first case, the temperature and specific humidity of ECMWF were linearly interpolated, while pressure was logarithmically interpolated to the GPS height.For the latter case, the temperature was extrapolated using the mean temperature gradient of the three lowest ECMWF layers.The pressure was calculated by stepwise application of the barometric height formula for each 20 m, while the specific humidity is estimated in parallel assuming the mean relative humidity from the two lowest ECMWF levels.

Summary of the uncertainty of the GNSS-derived IWV
We now can calculate the total uncertainty of the GNSSderived IWV after substituting Eqs. ( 25) and ( 27) to Eq. ( 16): Table 4 summarises the calculated total uncertainties of the GNSS-derived IWV for three GRUAN sites: LDBO, LDRZ, and NYA2.For each site, the GPS data acquired from the year 2014 were processed using a PPP strategy to obtain ZTD time series.The corresponding total ZTD uncertainties were then determined using the approach discussed in Sect.3.1.6and the mean values over the year were used.The estimated ZTD was converted to the IWV using the measured ground pressure and the mean temperature obtained from the ECMWF reanalysis data, ERA-Interim (Dee et al., 2011).The values of the conversion factor Q were calculated using Eq. ( 26).In Table 4, the corresponding absolute values for IWV, ZTD, ground pressure, and mean temperature are given using the mean values of the year 2014 for each site.As shown in Table 4, the uncertainties in the ZTD dominate the error budget of the resulting IWV, contributing over 75 % of the total IWV uncertainty.The impact of the uncertainty associated with the conversion factor Q increases slightly when the weather conditions get moist.
The IWV uncertainty is calculated for each data point in the GRUAN data products.Therefore, we first calculated the time series of the total ZTD uncertainty using the method discussed in Sect.3.1.1.Thereafter, the uncertainties for the ZHD and the conversion factor Q were determined using the uncertainties listed in Table 4 for each input variable.Finally, the total IWV uncertainty was calculated using Eq. ( 29).One example of the calculated total IWV uncertainties is shown in Fig. 9 for the month of June 2014.The mean IWV uncertainties are 0.68, 0.65, and 0.53 kg m −2 for LDB0, LDRZ, and NYA2, respectively.The corresponding variations of IWV uncertainties, given by the standard deviations, are 0.16, 0.13, and 0.06 kg m −2 .As expected, the variation in the IWV uncertainties is highly correlated to the variation in the ZTD uncertainties, and the larger uncertainties, as discussed in Sect.3.1.1,are caused both by the larger radial orbit errors and by a worse geometry of the satellites (less number of satellites visible).a The values are given by the mean ZTD uncertainty calculated from 1 year of data for LDB0, LDRZ, and NYA2, respectively.b For GRUAN sites equipped with surface barometers which are calibrated routinely.c Taken from Wang et al. (2005) based on the comparison between ECMWF reanalysis and radiosonde data.d Taken from Table 1 in Bevis et al. (1994).e Percentage of the total IWV uncertainty.f The constant given in Eq. ( 23).
As discussed above, currently the total IWV uncertainty obtained from a theoretical analysis is calculated by ignoring the site-dependent effects, i.e. signal multipath.If the GRUAN site is co-located with other techniques, the total IWV uncertainty can be estimated from the statistical method.One example is seen for the IGS station ONSA at the Onsala Space Observatory.The mean total uncertainty obtained from the theoretical analysis is 0.59 kg m −2 , while the uncertainty, given by a statistical analysis, presented in Table 1 is 0.7 kg m −2 when assuming no bias in the VLBI data.We interpret the difference to be mainly due to sitedependent effects ignored in the theoretical analysis.

Conclusions
Two methods were discussed in order to determine the total uncertainty of the GNSS-derived IWV.When there are at least three co-located techniques available, measuring the variability of the IWV at the same time, a statistical analysis is applied.This method is, however, difficult to apply on the present observational network because three independent methods for IWV measurement are not available.Therefore, a theoretical analysis, where the total uncertainty of the IWV is calculated from each one of the input variables according to the rule of uncertainty propagation for uncorrelated errors, is used.
In order to minimise the ZTD uncertainty caused by the factors discussed in Sects.3.1.2-3.1.5,some conditions (see Sect. 3.1.6)need to be fulfilled in the GPS data processing.The method discussed in Sect.3.1.1was applied to calculate the ZTD errors caused by the orbit errors for each time epoch which were added to the corresponding formal error.The calculated total ZTD uncertainties for three GRUAN sites are at the same level as the 1σ ZTD uncertainty (4 mm) calculated by the IGS.The uncertainty of the total IWV is dominated by the uncertainties in the ZTD, which contributes more than 75 % of the total IWV uncertainty.
The theoretical method will be implemented in the GRUAN GNSS central data processing.In summary, the following steps will be taken to calculate IWV uncertainty for each data point.Firstly, the ZTD uncertainty, including the systematic satellite orbit error and the formal error, is calculated.Secondly, the ZHD uncertainty is obtained using the uncertainty of ground pressure, estimated using the method presented by Immler et al. (2010), and the uncertainty in Eq. ( 23).Thirdly, the uncertainty for the conversion factor Q is calculated using the uncertainties of the mean temperature, k 3 , and k 2 listed in Table 4. Finally, the total IWV uncertainty for each data point is calculated using the uncertainties for the ZTD, the ZHD, and the conversion factor Q using Eq. ( 29).
For sites where two additional independent techniques are available, the IWV uncertainty estimated from the statistical method can be used to assess the stability of the data quality and potentially improve the operational theoretical method.For example, the statistical method can be used to quantify the ZTD uncertainty, which may change due to sitedependent effects such as signal multipath.
Although the method presented in the work is based on a PPP analysis, it can be applied to calculate uncertainties of the ZTD obtained from a network solution, but with some modifications.In PPP, the receiver clock parameters can partly compensate for the orbit errors.It is, however, not the case for a network solution where both satellite and receiver clock errors are cancelled out.Meanwhile, the baseline length and orientation in a network solution have a significant impact on the resultant ZTD uncertainties.

Figure 1 .
Figure 1.The impact of the radial and the tangential orbit errors in the estimated ZTD.

Figure 2 .
Figure 2. The ZTD error due to orbit errors for three GRUAN sites: (a) LDB0, (b) LDRZ, and (c) NYA2, where the blue dots show the simulated ZTD error, the red circles show the theoretical ZTD error caused only by the radial orbit error, and the green squares show the theoretical ZTD error caused only by the tangential orbit error.

Figure 3 .
Figure3.The impact of the elevation cutoff angle on the estimated ZTD for five GRUAN sites.The result obtained from the 5 • elevation cutoff solution is used as the reference, which however does not necessarily mean that it is the most accurate.Ideally, if all models are correct, the result shall not depend on the cutoff angle at all.

Figure 5 .
Figure 5.The time series of the estimated ZTD and the corresponding total ZTD uncertainty, for one day (1 June 2014) and for the three GRUAN sites: LDB0, LDRZ, and NYA2.

Figure 6 .
Figure 6.Values of f (λ, H ) calculated from different latitudes and from four different site heights.

Figure 7 .
Figure 7. Differences in the surface pressure measured by three different barometers.

Figure 8 .
Figure 8.The uncertainty of the conversion factor Q as a function of the uncertainty in the mean temperature (T m ) and for three groups of the uncertainties in the constants k 3 and k 2 : the nominal value (blue line), double of the nominal value (red line), and half of the nominal value (green line) for (a) a fixed uncertainty of k 2 and (b) a fixed uncertainty of k 3 .

Figure 9 .
Figure 9.The estimated total IWV uncertainty for the month of June 2014 and for three GRUAN sites: (a) LDB0, (b) LDRZ, and (c) NYA2.
The uncertainty of the ZWD only addressing the random errors.c The bias in the ZWD estimates for each technique.d The total uncertainty of the ZWD addressing both the random and systematic errors.e a The corresponding values were taken from Table2(synchronisation to all data) inNing et al. (2012).b
into account and assuming that all orbit errors are compensated by the ZTD.The other was calculated in the same way but only taking the tangential orbit errors into account.The results indicate that most portion of the orbit errors (> 90 %) are compensated by other two parameters (receiver clock and

Table 3 .
Mean and standard deviation (σ ) of the slant delay errors * .Corresponding values for M 0 , σ 0 , M , and σ used in Eqs. ( a

Table 4 .
Uncertainties in the GNSS-derived IWV calculated from the uncertainties associated with input variables.