Development of an improved empirical model for slant delays in the troposphere (GPT2w)

Global pressure and temperature 2 wet (GPT2w) is an empirical troposphere delay model providing the mean values plus annual and semiannual amplitudes of pressure, temperature and its lapse rate, water vapor pressure and its decrease factor, weighted mean temperature, as well as hydrostatic and wet mapping function coefficients of the Vienna mapping function 1. All climatological parameters have been derived consistently from monthly mean pressure level data of ERA-Interim fields (European Centre for Medium-Range Weather Forecasts Re-Analysis) with a horizontal resolution of 1°, and the model is suitable to calculate slant hydrostatic and wet delays down to 3° elevation at sites in the vicinity of the earth surface using the date and approximate station coordinates as input. The wet delay estimation builds upon gridded values of the water vapor pressure, the weighted mean temperature, and the water vapor decrease factor, with the latter being tuned to ray-traced zenith wet delays. Comparisons with zenith delays at 341 globally distributed global navigation satellite systems stations show that the mean bias over all stations is below 1 mm and the mean standard deviation is about 3.6 cm. The GPT2w model with the gridded input file is provided at http://ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/GPT2w/.


Introduction
Troposphere delay modeling is one of the main error sources in the analysis of space geodetic techniques operating at microwave frequencies, such as global navigation satellite systems (GNSS), very long baseline interferometry (VLBI), or Doppler orbitography and radiopositioning integrated by satellite (DORIS); see Nilsson et al. (2013) for a detailed overview. Normally, the concept of mapping functions is used to account for troposphere delays in the form of DL ðeÞ ¼ DL z h Á mf h ðeÞ þ DL z w Á mf w ðeÞ ð 1Þ The total slant delay DL at an elevation angle e is the sum of a hydrostatic and a wet portion, and each of them can be expressed as the product of a zenith delay and the corresponding mapping function. Hydrostatic zenith delay values DL h z are determined with sufficient accuracy from the local instantaneous pressure and the approximate station coordinates following Davis et al. (1985). In contrast, characteristic ''wet'' surface parameters, such as the water vapor pressure, represent largely inadequate measures of the non-hydrostatic vertical refractivity as required in evaluating zenith wet delays DL w z . These values are therefore usually estimated in the analysis of space geodetic observations as unknown parameters.
However, several positioning and navigation tasks like real-time applications do not have the benefit of post-processing analyses, necessitating the availability of accurate a priori estimates for DL w z . Widely adopted utilities including proxies for the zenith wet delays are RTCA-MOPS (1999), originally called UNB3 (Collins et al. 1996), the European Space Agency Galileo User Receiver model (ESA GAL-TROPO) (Krueger et al. 2004(Krueger et al. , 2005Martellucci 2012) and the concurrently developed TropGrid model (Krueger et al. 2004) with a recent update called TropGrid2 by Schüler (2014). The ESA model and both TropGrid versions account for annual and diurnal variations of the underlying parameters. Here, we present a new empirical zenith wet delay model, which, in combination with a corresponding utility for the pressure and fully consistent mapping functions, provides troposphere delays down to elevation angles of 3°after being fed with positional information and the time of observation.
A wide range of mapping functions has been used in the past. For example, the new mapping functions (NMF ;Niell 1996) and the global mapping functions (GMF; Böhm et al. 2006a) are empirical or the so-called blind mapping functions which only need the day of year and approximate station coordinates as input. Alongside such climatological approaches, the isobaric mapping functions (IMF; Niell 2001) and the Vienna mapping functions 1 (VMF1; Böhm et al. 2006b) account for the actual refractivity being derived from the operational analysis fields of numerical weather models at the epoch of the observations. Consequently, VMF1 positively exceeds NMF and GMF in terms of accuracy, but it is also prone to being inapplicable if the underlying mapping functions have not been updated.
Pressure values as required for the determination of zenith hydrostatic delays should be preferably taken from local measurements close to the antenna or from the gridded output of numerical weather models as, e.g., provided with the VMF1 at http://ggosatm.hg.tuwien.ac.at/. In case these streams are unavailable, analysts need to resort to blind models by Berg (1948) and Hopfield (1969), to UNB3m (Leandro et al. 2006), or to the global pressure and temperature model (GPT; Böhm et al. 2007) for approximate surface pressure information. Lagler et al. (2013) improved and combined GPT and GMF, then calling the updated blind model GPT2. They used 10 years (2001-2010) of 37 monthly mean pressure level data from the ECMWF (European Centre for Medium-Range Weather Forecasts) Re-Analysis (ERA-Interim; Dee et al. 2011) to determine mean values (A 0 ) as well as annual (A 1 , B 1 ) and semiannual amplitudes (A 2 , B 2 ) for selected parameters r on a regular 5°grid at mean ETOPO5 (earth topography) heights following where doy is the day of the year. GPT2 (Lagler et al. 2013) provides blind values of the hydrostatic and wet mapping function coefficients a h and a w , pressure p, temperature T and its lapse rate dT, and the water vapor pressure e. The underlying routine evaluates (2) at the four grid points surrounding the target location before extrapolating the parameters vertically to the desired height and interpolating the data from those base points to the observational site in horizontal direction. It should be noted here that the extrapolation of a h (strictly speaking of the hydrostatic mapping function) follows Niell (1996), whereas a w is assumed to be constant in the vicinity of the earth surface. The extrapolation of the pressure relies on an exponential trend coefficient related to the inverse of the virtual temperature (see Böhm et al. 2013, Eq. 23-27, ibid.), and the linear extrapolation of the temperature utilizes the GPT2 inherent temperature lapse rate dT. Surface grids for specific humidity within that model have been derived from linear interpolation between pressure levels in the vicinity of earth's surface. These parameters could be possibly used to determine values of zenith wet delays, e.g., by using the expressions of Saastamoinen (1972), though this approach is not optimal and in the following we are introducing GPT2w as an extension to GPT2 with an improved capability to determine zenith wet delays in blind mode. The next describes the development of this new, auxiliary model, which is then validated against zenith delays from GNSS observations. An outlook addresses possible applications of GPT2w and concludes the study.

Development of GPT2w
The zenith wet delay expressions as provided by Askne and Nordius (1987;Eq. 18, ibid.) represent the starting point of our derivations, Herein, k 2 0 and k 3 are empirically determined coefficients, R d denotes the specific gas constant for the dry constituents, g m is the gravity acceleration at the mass center of the vertical column of the atmosphere, and e s stands for the water vapor pressure at the site. The weighted mean temperature at height H along the local vertical dz can be formulated as and the water vapor decrease factor k is defined implicitly via the expression using the relationship between surface pressure values (e s plus the total surface pressure p s ) and those defined on any target level (e, p). In essence, to use (3) for the determination of zenith wet delays, we need the water vapor pressure at the site, the weighted mean temperature T m , and the water vapor decrease factor k. While T m is readily accessible via numerical integration through the pressure level data, the derivation of k evades a straightforward treatment because it can be calculated for any pair of pressure levels, always leading to significantly different results which is due to the variability of water vapor with height. On the other hand, we are able to determine the zenith wet delays by numerical integration of the wet refractivity along the site vertical ), allowing us to simply invert (3) toward a global k grid, which will be very much representative of the decrease factor behavior through the entire troposphere and fully consistent with the zenith wet delay. Figure 1 illustrates mean values of k, its annual and semiannual amplitudes, and the standard deviations of the residuals as estimated by least-squares adjustment of the ''observed'' water vapor decrease factor field obtained by inverting (3). Gauged to 10 years of monthly mean pressure level data from ERA-Interim (1°horizontal resolution) via DL w z and T m , the fitted k grids are fully consistent with GPT2. Both the global mean field (Fig. 1a) and the annual cycle ( Fig. 2a) show typical atmospheric circulation patterns, e.g., the intertropical convergence zone standing out as a belt of small vertical gradients in the wake of upwelling moist air (Wallace and Hobbs 1977). On the contrary, horse latitudes (at about 30°north and south) are characterized by a large-scale subsidence of dry air, which in combination with strong regional evaporation over oceanic surfaces (e.g., Arabic and Mediterranean Sea, west coasts of Africa and America) implicates steep gradients of k. Overall, Fig. 1 clearly demonstrates that it is not sufficient to apply constant decrease factors for the water vapor pressure, neither in space nor in time.
Analogous plots of mean values, annual and semiannual amplitudes, as well as post-fit standard deviations are depicted in Fig. 2 for the weighted mean temperature T m of the water vapor pressure (4). The distribution of the mean values is mostly latitude-dependent, while a dominant annual variation of 10-15 K prevails over north-east Asia and the eastern part of Canada. Semiannual amplitudes of T m are strongest at latitudes higher than 60°and in subtropical regions experiencing bimodal (two-peak per year) patterns in precipitation and temperature attached to two distinct wet seasons; see, e.g., northern India, the Arabic Sea, and the western Sahel zone. The implications of the above-described spatial patterns shall be exemplified for one particular grid point, chosen to be located in the Arabic Sea (latitude 22.5°north, longitude 63.5°east), where a prominent water vapor decrease factor has been noted (Fig. 1a). In Fig. 3, we plot the values of water vapor pressure, weighted mean temperature, water vapor decrease factor, and the zenith wet delay as derived from the monthly mean ERA-Interim fields and as approximated with GPT2w over 3 years (2001.0-2003.12). The fitted climatological mean in GPT2w gives a particularly good account of the actual annual variation over the selected time span, while larger residuals remain for the semiannual cycle; see, e.g., Fig. 3c and recall the standard deviations displayed in Fig. 1d. Interestingly, the broader annual maxima in the local water vapor pressure values and the strong semiannual pattern of k act to generate distinct and sharp annual peaks in the zenith wet delay of about 10 cm.
Having shown the relevance of a precise consideration of T m and k for zenith wet delays, we add these two prime parameters to our set of variables determined from 10 years of ERA-Interim data. The horizontal resolution of 1°in latitude and longitude is not reduced toward coarser spacing as in the case of GPT2 because the wet part in the atmosphere has more small-scale structures, in particular at coastal areas (Fig. 1). With GPT2w, the vertical extrapolation of e s adheres to (5), utilizing the GPT2w-inherent values of the water vapor decrease factor k. Pressure reductions conform to GPT2's use of exponential trend coefficients calculated from grid point-wise virtual temperature information. Such isothermal scale heights may be alternatively adjusted for adiabatic effects, but the benefit of this approach on the side of delay predictions is not entirely conclusive (Schüler, 2014). Table 1 summarizes the main features of the newly suggested model, which can be compared to Tab. 1 by Lagler et al. (2013).

Validation with zenith total delays from GNSS
Since 2011, the United States Naval Observatory (USNO) archives and distributes final tropospheric estimates from the observation data of all GNSS sites of the International GNSS Service (IGS) in the Solution Independent EXchange format (SINEX) with a temporal resolution of 5 min and with a latency of about 4 weeks (Hackman and Byram 2012). The accuracy of the tropospheric delays in zenith direction is specified with 4 mm by the IGS Central Bureau at http://igs.org/components/prods.html. For the purpose of validating GPT2w, all available SINEX files in 2012 (about 110,000) were downloaded and subsequently cleansed. A visual inspection showed that in some cases, no realistic tropospheric estimates could be derived from the GNSS observations and we consequently excluded these specific records from the comparison in order to avoid falsification of the validation outcome. Furthermore, only zenith total delays (i.e., the sum of hydrostatic and wet zenith delays) with a formal error smaller than 18 mm were retained, which is a reasonable threshold yielding a total of 341 GNSS sites with at least 110 days of tropospheric delays. In total, we removed about 3.4 % of all entries of zenith total delays in the SINEX files.  Using the positional information and the (mostly 6 hourly) temporal sampling of all 341 GNSS sites, we deployed GPT2w to obtain blind predictions of zenith total delays. The approach by Davis et al. (1985) was applied to determine the zenith hydrostatic delays, and we evaluated (3) to estimate the zenith wet delays. Figure 4 illustrates station-wise mean values (i.e., biases) and empirical standard deviations (i.e., RMS, root-mean-square quantities after removing the biases) of the differences between the IGS and GPT2w zenith total delays series.
The overall agreement between the GPT2w-based zenith total delays and those from IGS observations is in the range of a few centimeters. Biases (Fig. 4, top) vary between -4.2 cm and ?7.3 cm, exhibiting small but persistently negative values for Europe but generally positive values over North America. Largest offsets occur for small islands (e.g., Hawaii or station DGAR, Diego Garcia Island, central Indian Ocean), where the horizontal resolution of 1°is apparently not sufficient to capture microclimate and where the station heights do not agree well with the mean ETOPO5 topography. Sites with short record lengths down to 30 % of the entire year 2012 (in the Caribbean and Middle America) quite likely hold unrepresentative bias estimates but are nonetheless included in Fig. 4a. RMS differences (Fig. 4, bottom) with IGS data cover a range of 1.5 to 6.1 cm and are of typical low magnitude for dry (sub-) polar regions, innercontinental areas, and high-terrain stations such as the two North-Andean receivers located above 2,500 m in Colombia and Ecuador. Significant non-climatological variability in zenith total delays, presumably related to westerly cyclone-anticyclone patterns or monsoon activity, produces large RMS values beyond 4 cm for 100 midlatitude and subtropical stations at an average height of 149 m. The confinement of such peak disparities to altitudes near sea level has already been exemplified by akin comparisons in Schüler (2014). Despite the obvious limitations of our blind model, the majority of sites in the Asian-Pacific region profit from using GPT2w instead of GPT2 (with low-spatial resolution and inadequate wet delays) or the ESA model that disregards semiannual harmonics. Prime examples where the inclusion of semiannual terms allows for a much better representation of non-sinusoidal delay variability over 1 year are stations PBR2 (Port Blair, Bay of Bengal) and LCKI (Lucknow, Indian mainland); see Fig. 5.
To further substantiate the enhanced performance of GPT2w in predicting the IGS zenith total delays with respect to other models, we determined bias plus RMS statistics at the 341 test locations for RTCA-MOPS (1999), the ESA blind model (ESA GALTROPO; Krueger et al. 2004Krueger et al. , 2005Martellucci 2012), and GPT2 (Lagler et al. 2013) by analogy to Fig. 4. The global averages of these values for each model are summarized in Table 2. Marked biases in RTCA-MOPS (-25 mm) and ESA (8 mm) have been reduced to only a slight overestimation of the total delay by GPT2 (-3 mm), a residual that is excellently accounted for by the improved zenith wet delay formulation within GPT2w. RMS values are at a level of 3.8 cm both for the ESA model and GPT2, but drop to 3.6 cm when utilizing GPT2w. As surmised above, the distinguished ability of GPT2w to represent the temporal variability in the observations mainly resides in its use of semiannual harmonics, which accounts for 1.6 mm of the RMS reduction in Table 2.
Of course, a 1-year analysis comparing predicted zenith delays against those from the IGS is rather tentative and should be extended in future. Schüler (2014; Fig. 1, ibid.) provides statistics for several years illustrating that the RMS values change by about 2 mm over the years depending on the distribution of IGS antennas available. Considering this distribution of RMS values, our results agree reasonably well with those for TropGrid2 by Schüler (2014). The RMS values in his study, however, contain systematic contributions from the biases between the time series. A telling graphical illustration of individually improved statistics is presented in Fig. 6 in the form of station-wise differences in RMS values with respect to the IGS zenith total delays. Residuals reckoned in the sense ''GPT2w minus ESA'' and ''GPT2w minus GPT2'' are almost continuously of negative sign (288 instances for ESA, 253 instances for GPT2), suggesting that GPT2w supplies more accurate approximations of the observed delays. Maximum reductions greater than 2 cm are apparent for the two receivers operated at LCKI and LCK2 (both at Lucknow, India). GPT2w acts to redress specific shortcomings of GPT2 in equatorial regions and the southern subtropics, whereas advances with regard to the ESA model are particularly evident for northern hemispheric midlatitudes.
An extended validation of GPT2w with IGS data and ray-traced delays as well as a comparison with state-of-theart tropospheric models like RTCA-MOPS or the ESA tropospheric blind model can be found in Möller et al. (2014). The bias and RMS statistics in that research were, however, computed from the globally merged sample of zenith total delay differences instead of determining them as mean values over all stations as above. Furthermore, the grid files of water vapor pressure have been slightly modified in the course of the present study.

Summary and outlook
With GPT2w, we have introduced a new blind tropospheric delay model which is based on gridded values of water vapor pressure, water vapor decrease factor, and weighted mean temperature. In terms of zenith total delays, the globally averaged bias is below 1 mm and the RMS difference is about 3.6 cm as when compared to zenith total delays from GNSS at 341 globally distributed sites. Since GPT2w is also equipped with fully consistent hydrostatic and wet VMF1 coefficients, it may not only be used for positioning and navigation purposes but also for highprecision applications, like geophysical studies, where the  wet mapping functions are essential to estimate residual zenith wet delays. In the future, the addition of diurnal and semidiurnal amplitudes will be considered for selected parameters. GPT2w also contains the mean values as well as annual and semiannual amplitudes of the weighted mean temperature. This is an important quantity for the determination of the integrated water vapor or precipitable water as required in GNSS meteorology (Bevis et al. 1992). However, detailed studies on that application still have to be carried out.
Another promising future application of GPT2w is the combination with local meteorological observations at the sites. It will be investigated whether, e.g., a local measurement of water vapor pressure can be reasonably connected to climatological values of weighted mean temperature and water vapor decrease factor in order to determine improved zenith wet delays or whether additional corrections will have to be applied to fully exploit local measurements.
The MATLAB source code of GPT2w is available at: http://ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/GPT2w/. Johannes Böhm is Professor for Space Geodetic Techniques at the Department of Geodesy and Geoinformation at the Vienna University of Technology, Austria. In his PhD studies and thereafter, he dealt with tropospheric path delays and developed models for the analysis of observations from space geodetic techniques. His main interests now are atmospheric effects in space geodesy in general and the analysis of very long baseline interferometry observations. Gregor Möller is University Assistant and PhD student at Vienna University of Technology. His current research focuses on GNSS signal analysis and troposphere delay modeling in terms of meteorology and navigation applications. Since 2011, Gregor Möller is a member of the IGS Troposphere Working Group. For his work in the field of troposphere delay modeling, he received ION's Best Presentation Award in September 2013.
Michael Schindelegger is University Assistant at Vienna University of Technology and recently completed his PhD studies in the field of Earth rotation and its excitation from geophysical fluids. His core interests are atmospheric circulation and ocean tides, troposphere delay modeling, and the usability and reliability of global analysis models for space geodetic applications.
Gregory Pain is a graduate of the IGN engineering school. During a traineeship in summer 2013, he visited the Vienna University of Technology.

Robert Weber is Associate
Professor at the Department of Geodesy and Geoinformation, TU Vienna, and former chair of the IGS GNSS Working Group. His main fields of research include geodetic reference systems, active GNSS reference station networks, and the applications of GNSS for geodynamics and meteorology.