Tropospheric delay performance for GNSS integrated water vapor estimation by using GPT2w model, ECMWF's IFS operational model and in situ meteorological data

2.1 GNSS and meteorological datasets

In this study, GNSS data from two permanent stations from Greece and Cyprus were used. The first is the Elassona station, named ELAS, located at Elassona in Thessaly, Central Greece (39.53^∘ N, 22.12^∘ E, 303 m) managed by the GNSS_QC team of the Aristotle University of Thessaloniki (Fotiou and Pikridas, 2012; Pikridas et al., 2014). The second is the Nicosia station located at Athalassa to the southeast of Nicosia, named NICO (35.08^∘ N, 33.23^∘ E, 162 m) that belongs to the Cyprus Department of Meteorology; this station participates in the EPN/EUREF network (Bruyninx, 2004). Elassona station is equipped with Leica receiver (model GR30) and is collocated with Vaisala meteostation (model PTU300), whilst the Nicosia station is also equipped with Leica receiver (model GR25) and is located 15 m from the main synoptic meteorological station of Athalassa (WMO17607) with a proprietary data acquisition system utilising Vaisala meteosensors. The temperature, relative humidity and pressure sensors provide measurements every 2 s and the values are recorded as a 10 min linearly averaged value (World Meteorological Organization, 2010). GNSS stations have an observation recording rate of 30 s and the capability to record all available satellite constellations' signals (GPS, GLONASS, GALILEO, BEIDOU). Figure 1 shows the cluster definition of permanent GNSS network including IGS and EPN/EUREF stations that was employed during data processing with Bernese software.

A time span of 8 days was selected, namely, from 27 May until 3 June 2018. Under the frame of the European research project BeRTISS (Balkan-Mediterranean Real Time Severe weather Service), the GNSS data were analyzed on an hourly basis. The data were recorded using a satellite elevation cut-off angle of 10^∘. It is important to note that in order to derive absolute ZTD estimates, the average length of baselines included in the processing should be over 500 km which coincides with the selection of several EUREF and IGS stations. These stations are included in the network processing scheme for datum definition and consistent absolute tropospheric estimation. Long baselines reduce the high level of correlation of the tropospheric delays which occurs on short baselines (Duan et al., 1996). The processing scheme concerns the use of the Bernese software v5.2 developed by the Astronomical Institute of the University of Bern (Dach et al., 2015). The ultra-rapid precise orbit information which refers to the current reference frame is used. IGS has decades of products that are archived which can be accessed through the ftp. To reduce the age of the prior, discontinued Predicted orbits, the IGS started the Ultra-rapid products officially in November 2000 (GPS week 1087). Ultra-rapid products are available for real time and near real time use. The products are released four times per day. The IGS_08.atx model with absolute antenna calibration values was applied. The Saastamoinen model (see Saastamoinen, 1972) with the Global Mapping Function (GMF) mapping function (Böhm et al., 2006) was also used in the processing. Once of the coordinates of the fixed station were set, the coordinates for the rest of the stations were obtain by means of the maximum observations percentage between the fixed stations and the others. Solutions for each processing sessions were combined with ADDNEQ algorithm of Bernese software and the final solution was derived. All the initial phase ambiguities of carrier frequencies were resolved using the QIF (Quasi Ionosphere Free) strategy (Dach et al., 2015). As it is known, carrier phase data are biased by an integer number of carrier wavelengths which called initial phase ambiguity and must be estimated from the data. In case phase ambiguity resolution is succeeded, high accuracy results are obtained. As a consequence, ZTD estimates were obtained with an accuracy of a few (2–3) millimetres.

Figure 1Map of GNSS network including IGS and EPN/EUREF stations that was employed during data processing.

Download

2.2 Retrieval of IWV from GNSS ground receivers – Converting the estimated zenith wet delays onto precipitable water

The delay which the wet component of the troposphere induces on GPS signals provides also an opportunity for sensing Precipitable Water (PW) with ground-based stations (Bevis et al., 1994). Precipitable Water (PW) is defined as the depth of water that would result if all atmospheric water vapor in a vertical column of air is being condensed to liquid. For this reason, Zenith Wet Delays (ZWD) estimated by GNSS can be converted to PW using the formula (Bevis et al., 1994):

where, the ZWD is derived from the relation $\text{ZWD}=\text{ZTD}-\text{ZHD}$, with ZHD depending on the dry air gases accounting for the greatest part of the tropospheric delay and

where, $k_{\mathrm{2}}^{\prime }=k_{\mathrm{2}}-m{k}_{\mathrm{1}}$, k₁=77.604 K hPa⁻¹, k₂=17 K hPa⁻¹ and $k_{\mathrm{3}}=\mathrm{3.776}\times {\mathrm{10}}^{\mathrm{5}}$ K hPa⁻¹ and ρ=1025 kg m⁻³ is the density of liquid water, R_u=461.51 $\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{kg}}^{-\mathrm{1}}$ is the specific gas constant of water vapor, T_m is the water vapor weighted mean temperature of the atmosphere as defined by Davis et al. (1985) and m is the ratio of the molar mass of water vapor and dry air. The physical constants k₁, k₂ and k₃ are from the well-known formula for atmospheric refractivity of Davis et al. (1985). In general, PW is related to ZWD by a factor (Π), that is approximately 0.15, however the actual value of factor Π can vary by as much as by 20 % (Bevis et al., 1994; Rózsa, 2014) since it is a a function of the weighted mean temperature of the atmosphere which varies with location, altitude, season and weather.

As a result, in calculating PW from Eqs. (1) and (2), the largest source of error is attributed to the mean temperature which varies according to location, height, weather and season. More specifically, the uncertainties in Π derive from the uncertainties in the mean temperature of the atmosphere T_m and in the physical constants k₁, k₂ and k₃. Therefore, factor Π is a function of the weighted mean temperature of the atmosphere (Davis et al., 1985) and can be determined to about 2 % when it is computed as a function of surface temperature, or 1 % if data from numerical weather models are used. (Bevis et al., 1994; Businger et al., 1996). As it concerns the parameter T_m, Bevis et al. (1992) suggested that it can be estimated from the surface temperature. In our study, the mean temperature of the atmosphere in Kelvin, T_m, was estimated according to the formula by Mendes et al. (2000) which holds for mid-latitudes:

where, T_s is the surface temperature in Kelvin.

In this analysis, we adopted two ways for estimating the PW values (conversion to IWV). First, we retrieved (after routine analysis) the associated parameters from GPT2w model in order to compute the tropospheric component ZHD. Based on the approximation by Saastamoinen (1972), the ZHD can also be expressed as:

where, P is the atmospheric pressure, Φ is the station latitude and h the station's orthometric height. Subsequently, we compute the ZWD and finally the quantity IWV (Pikridas et al., 2014). The second calculation scenario was based on the estimation of the same parameters but instead of GPT2w values, we used the ground (real) meteorological records obtained from Vaisala sensors collocated to the respective GNSS receivers at Nicosia and Elassona sites.

2.3 GPT2w (Global pressure and temperature 2 wet) model data

The blind empirical model GPT2w (Böhm et al., 2014) produces the mean value plus annual and semiannual amplitudes of the pressure among else. Firstly, Böhm et al. (2007) constructed a blind model, named GPT, based on 3-year monthly mean profiles of pressure and temperature from the ECMWF 40-year reanalysis data (ERA40) for giving pressure and temperature for any user position and epoch (i.e., the measurement interval of a GPS receiver). Subsequently, Lagler et al. (2013) created GPT2, by improving and combining GPT and Global Mapping Function (GMF) using 10 years (2001–2010) of 37 monthly mean pressure isobaric level data from the ECMWF reanalysis (ERA-Interim). The latest GPT2w is an upgrade of GPT2 with an improved capability in the calculation of zenith wet delays in the blind mode. The GPT2w model with the gridded input file is available at http://ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/, last access: June 2018. GPT2w and is based on global $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ gridded values of tropospheric variables such as surface pressure for modelling the zenith hydrostatic delay or water vapor pressure, weighted mean temperature and water vapor decrease factor for modelling the zenith wet delay. To determine the pressure and temperature values at a GNSS station, the model first calculates the pressure and temperature values at four nearest grid points surrounding the station. Then, temperature is corrected for the altitude difference between the station and grid. For this reason, at each of the four grid points, the altitude gradient and the temperature lapse rate are added to the temperature value. As regards the pressure, values at these four grids at the height of the station are calculated using a vertical exponential trend coefficient related to the inverse of the virtual temperature. Finally, the pressure and temperature at the GNSS station are interpolated from the above four respective grid values. Nevertheless, according to several studies, GPT2w is considered a very good blind model, in general (Liu et al., 2017).

2.4 ECMWF's Integrated Forecasting System (IFS) operational model datasets

In this study, Zenith hydrostatic delays (ZHD) and non-hydrostatic (wet) delays (ZWD) are also obtained from the Potsdam mapping functions (PMFs) (Zus et al., 2014; Balidakis et al., 2018) with ECMWF's Integrated Forecasting System Cycle 43r1 featuring as the underlying numerical weather model. For the ray-tracing, the refraction index tensor and the position operator thereof are composed of 6-hourly employing surface pressure and geopotential fields at a horizontal resolution of approximately 9 km (O1280 spectral), and temperature and specific humidity tensors discretized at 137 model levels (L137). The IWV is calculated by numerical integration in the model level of ECMWF operational analysis, starting from the point of interest, which is the reference point of the GNSS station (NICO, ELAS), to the top of the model (∼80 km). Our choice of ECMWF's IFS is deliberate for it is a state-of-the-art model repeatedly tested for geodetic and remote sensing purposes. The data assimilation is based on a 12-hourly weak constraint four-dimensional variational system employing orography-following hybrid sigma-pressure levels. The high spatial resolution allows the description of orography-related fluctuations in the atmospheric state that in datasets such as ERA-Interim (80 km) would remain unresolved. This is particularly relevant for sites situated in regions with steep orographic gradients such as ELAS.

Figure 2Variations of T, P, ZHD, ZWD, and IWV for Nicosia and Elassona stations from collocated meteorological stations (blue dots), GPT2w (red boxes) and ECMWF (black boxes) spanning the period 27 May–3 June 2018.

Download

Figure 3Bias between in situ T, P observations and T, P values from GPT2w (blue dot) and ECMWF models (red box). Also, bias between ZHD, ZWD, and IWV values retrieved from GNSS observations by using T, P from collocated meteorological stations and the respective ZHD, ZWD, IWV values retrieved from GNSS observations by using T, P from GPT2w model (blue dots), the respective ZHD, ZWD, IWV values retrieved from ECMWF model (red boxes) for Nicosia and Elassona stations during 27 May–3 June 2018.

Download

2.5 Methodology

To investigate the tropospheric delay performance for GNSS integrated water vapor estimation, we employed: (a) GPT2w model data, (b) ground meteorological records from meteo-stations collocated with the GNSS receivers and (c) ECMWF's (IFS) model data for the two sites during the period from 27 May to 3 June 2018, when abrupt heavy midday rainfalls, almost on a daily basis occurred at the station in Cyprus and fair weather conditions prevailed at the Greek station. It is worth noting that while the average area precipitation for the entire island is 19.6 mm for May and 13.2 mm for June, the local thunderstorms of 27–31 May and 1–3 June, produced an area averaged precipitation of 19.5 and 13.2 mm, respectively.

First, the capability of the blind model GPT2w and of the ECMWF (IFS) numerical prediction model to reproduce surface pressure and temperature variables is validated against in situ meteorological measurements over Nicosia and Elassona sites, by assessing the absolute bias (A1) and relative percentage deviations (D1) as follows:

In addition, the Standard Deviation and RMSE (root-mean-square error) of the differences between observed and modelled values for each station are computed and the time series of observed and modelled pressure and temperature variables are provided as well. The number of observations (points) used for RMSE are N=192 and N=65 for Nicosia and Elassona stations respectively.

Secondly, the tropospheric products: zenith hydrostatic delay (ZHD), zenith wet delay (ZWD) and integrated water vapor (IWV) are assed at the two sites as follows:

1. They are retrieved from GNSS observations by using surface temperature and pressure values from GPT2w model;

2. They are retrieved from GNSS observations by using surface pressure and temperature values from in situ meteorological stations;

3. By using ECMWF (IFS) operational model data as described at Sect. 2.4.

The results from the above three approaches are intercompared for both stations and for all three parameters (ZHD, ZWD, IWV), by means of time series and by calculating the absolute (A2, 3) and the relative (D2, 3) differences:

1. between parameters retrieved from GNSS observations by using GPT2w model atmospheric variables (T, P) and the respective ones retrieved from GNSS observations by using T, P obtained from in situ meteorological stations;

2. between parameters retrieved by using ECMWF (IFS) model and the respective ones retrieved from GNSS observations by using T, P obtained from in situ meteorological stations.

The following equations for the calculation of these differences are applied:

where, X represents either of ZHD, ZWD and IWV.