Inter-technique validation of tropospheric slant total delays

Atmos. Meas. Tech., 10, 2183–2208, 2017 https://doi.org/10.5194/amt-10-2183-2017 © Author(s) 2017. This work is distributed under the Creative Commons Attribution 3.0 License. Inter-technique validation of tropospheric slant total delays Michal Kaˇ cmaˇ rík1, Jan Douša2, Galina Dick3, Florian Zus3, Hugues Brenot4, Gregor Möller5, Eric Pottiaux6, Jan Kapłon7, Paweł Hordyniec7, Pavel Václavovic2, and Laurent Morel8 1Institute of Geoinformatics, VŠB – Technical University of Ostrava, Ostrava, Czech Republic 2Geodetic Observatory Pecný, Research Institute of Geodesy, Topography and Cartography, Zdiby, Czech Republic 3GFZ German Research Centre for Geosciences, Potsdam, Germany 4Atmospheric Composition Department, Royal Belgian Institute for Space Aeronomy, Brussels, Belgium 5Department of Geodesy and Geoinformation, Vienna University of Technology, Vienna, Austria 6Royal Observatory of Belgium, Brussels, Belgium 7Institute of Geodesy and Geoinformatics, Wrocław University of Environmental and Life Sciences, Wrocław, Poland 8GeF Laboratory, ESGT – CNAM, Le Mans, France Correspondence to: Michal Kaˇ cmaˇ rík ([email protected]) Received: 9 November 2016 – Discussion started: 5 January 2017 Revised: 13 April 2017 – Accepted: 3 May 2017 – Published: 12 June 2017 Abstract. An extensive validation of line-of-sight tropospheric slant total delays (STD) from Global Navigation Satellite Systems (GNSS), ray tracing in numerical weather prediction model (NWM) fields and microwave water vapour radiometer (WVR) is presented. Ten GNSS reference stations, including collocated sites, and almost 2 months of data from 2013, including severe weather events were used for comparison. Seven institutions delivered their STDs based on GNSS observations processed using 5 software programs and 11 strategies enabling to compare rather different solutions and to assess the impact of several aspects of the processing strategy. STDs from NWM ray tracing came from three institutions using three different NWMs and ray-tracing software. Inter-techniques evaluations demonstrated a good mutual agreement of various GNSS STD solutions compared to NWM and WVR STDs. The mean bias among GNSS solutions not considering post-fit residuals in STDs was −0.6 mm for STDs scaled in the zenith direction and the mean standard deviation was 3.7 mm. Standard deviations of comparisons between GNSS and NWM ray-tracing solutions were typically 10 mm ±2 mm (scaled in the zenith direction), depending on the NWM model and the GNSS station. Comparing GNSS versus WVR STDs reached standard deviations of 12 mm ±2 mm also scaled in the zenith direction. Impacts of raw GNSS post-fit residuals and cleaned residuals on optimal reconstructing of GNSS STDs were evaluated at intertechnique comparison and for GNSS at collocated sites. The use of raw post-fit residuals is not generally recommended as they might contain strong systematic effects, as demonstrated in the case of station LDB0. Simplified STDs reconstructed only from estimated GNSS tropospheric parameters, i.e. without applying post-fit residuals, performed the best in all the comparisons; however, it obviously missed part of tropospheric signals due to non-linear temporal and spatial variations in the troposphere. Although the post-fit residuals cleaned of visible systematic errors generally showed a slightly worse performance, they contained significant tropospheric signal on top of the simplified model. They are thus recommended for the reconstruction of STDs, particularly during high variability in the troposphere. Cleaned residuals also showed a stable performance during ordinary days while containing promising information about the troposphere at low-elevation angles. 1 Introduction Tropospheric slant total delay (STD) represents the total delay that undergoes the GNSS radio signal due to the neutral atmosphere along the path from a satellite to a ground receiver antenna. This total delay can be separated into the hydrostatic part, caused by the dry atmospheric constituents, and the wet part caused specifically by water vapour. By quantifying the total delay, and by separating the hydrostatic Published by Copernicus Publications on behalf of the European Geosciences Union. 2184 M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays and wet parts, it is possible to retrieve the amount of water vapour in the atmosphere along the path followed by the GNSS signal. During the processing of GNSS observations only the total delay in the zenith direction (zenith total delay, ZTD) above the GNSS antenna can be estimated for each epoch or for a time interval. ZTDs from GNSS reference stations are operationally assimilated into numerical weather prediction models (NWMs) for almost a decade (Bennitt and Jupp, 2012; Mahfouf et al., 2015). In Europe, this activity is coordinated mainly in the framework of the EUMETNET EIG GNSS Water Vapour Programme (E-GVAP, 2005–2017, phases I– III, http://egvap.dmi.dk). Many recent studies demonstrated a positive impact of the ZTD or integrated water vapour (IWV) assimilation on precipitation weather forecasts, especially of the short-time ones (Vedel and Huang, 2004; Guerova et al., 2006; Shoji et al., 2009; Guerova et al., 2016). In contrast, continuous developments in NWM forecasting and nowcasting tools, as well as increasing needs for better predictions of severe weather events, stress the demand of high-quality humidity observations with high spatial and high temporal resolutions. While ZTDs provide information in zenith directions above GNSS stations, linear horizontal tropospheric gradients give information about the first-order spatial asymmetry around the station. Besides, slant tropospheric delays (STDs) can provide additional details about the horizontal asymmetry in the troposphere, more specifically in the directions from a receiver to all observed GNSS satellites. With the increasing number of GNSS systems and satellites, the atmosphere scanning will be more complete, hence gaining even more interest. Bauer et al. (2011) showed a positive impact of STD assimilation into the Mesoscale Model 5 (MM5) and Kawabata et al. (2013) demonstrated a significant advantage of assimilating STDs into a high-resolution model in the case of forecasting local heavy rainfall event against the scenario of assimilating ZTDs only. Also, Shoji et al. (2014) and Brenot et al. (2013) showed promising techniques for prediction of severe weather events using advanced GNSS tropospheric products such as horizontal gradients and STDs. The GNSS tomography technique aiming at the three-dimensional reconstruction of the water vapour field (Flores et al., 2001) uses STDs as input data as well. Obviously, the quality of the tomography depends on both the accuracy of the STDs (Bender et al., 2009) and the observation geometry (Bender et al., 2011). Validation of GNSS slant delays with independent measurements is not a new research topic. GNSS slant delays were validated against water vapour radiometer (WVR) measurements in Braun et al. (2001, 2002) and Gradinarsky (2002). First attempts to derive slant delays from NWM fields and to compare them with GNSS STDs were carried out by De Haan et al. (2002) and Ha et al. (2002). Additional effort to evaluate GNSS slant delays using WVR and NWM data was done at GFZ Potsdam over the last few years. Bender et al. (2008) showed an existing high correlation within the three sources (GPS, WVR, NWM) of slant wet delays (SWDs) and tried to quantify the effect of removing multipath from GPS post-fit residuals using a stacking method what was also done by Kaˇ cmaˇ rík et al. (2012). Deng et al. (2011) validated tropospheric slant path delays derived from singleand dual-frequency GPS receivers with NWM and WVR data. Shangguan et al. (2015) compared GPS versus WVR slant IWV values (SIWVs) using a 184day dataset. They also analysed the influence of the elevation angle setting and the meteorological parameters (used for the conversion to IWV) on the comparison results. More recently, a validation of multi-GNSS slant total delays retrieved in real time from GPS, GLONASS, Galileo and BeiDou was presented by Li et al. (2015a) using WVR and NWM as independent techniques for the assessment. Using multiple GNSS constellations brought a visible advantage, in terms of not only the number of available slants but also their higher accuracy and robustness. Nevertheless, most of the studies presented thus far were limited to only a single strategy for obtaining GNSS STDs and usually restricted to a limited set of stations and/or a relatively short time period. The main purpose of this study is an extensive comparison of various solutions from GNSS processing, NWM ray tracing and WVR measurements using one common dataset as well as a comparison of results from collocated stations. The GNSS solutions evaluated in this work used 5 different software programs and 11 strategies and exploited the GNSS4SWEC benchmark dataset (Douša et al., 2016). Then, the paper studies the impact of various approaches on STD estimates and aims to find the most suitable strategy for estimating the GNSS-based STDs. Section 2 briefly introduces the validation study dataset, and Sect. 3 describes the process of retrieving GNSS STDs including an overview of the different GNSS solutions. Section 4 provides a description of STDs generated from NWMs, and Sect. 5 summarizes WVR principals and WVRbased STD solutions. Section 6 introduces the methodology used in the validation of STDs, and Sects. 7 and 8 study the results achieved at single GNSS reference stations and at closely collocated stations, respectively. 2 Experiment description The presented work has been carried out in the context of the EU COST Action ES1206 “Advanced Global Navigation Satellite Systems tropospheric products for monitoring severe weather events and climate” (GNSS4SWEC; http: //www.cost.eu/COST_Actions/essem/ES1206, 2013–2017). Three mutually cooperating working groups (WG) have been established to cover the proposed topics: (1) WG1 for Advanced GNSS processing techniques, (2) WG2 for GNSS for severe weather monitoring, and (3) WG2 for GNSS for climate monitoring. Validation of STDs belongs mainly under WG1, which is oriented toward the development of new Atmos. Meas. Tech., 10, 2183–2208, 2017 www.atmos-meas-tech.net/10/2183/2017/ M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays 2185 advanced tropospheric products, among other topics. The idea of preparing a common benchmark dataset, which could serve efficiently for most planned activities, was designed in the beginning of the project; the data were collected, cleaned and documented, and reference products were generated and assessed (Douša et al., 2016). The selected geographical area is situated in central Europe (Austria, Germany, the Czech Republic, Poland) where severe weather events, including extensive floods on Danube, Moldau and Elbe rivers, occurred between May and June 2013. The benchmark dataset gathers observations from 430 GNSS reference stations, 610 meteorological synoptic stations, 21 radiosonde launching sites, 2 WVR, 2 meteorological radars and output fields from the ALADIN-CZ NWM over a period of 56 days. ZTDs and horizontal tropospheric gradients from the reference GNSS and NWM-derived tropospheric products were already evaluated, and all resulted in very good agreement (Douša et al., 2016). All STDs used in this paper were computed by exploiting the benchmark dataset. From the complete benchmark dataset, we selected a subset of 10 GNSS reference stations situated at six different locations (Table 1). The selection was based on the following requirements: (1) long-term quality of observations and its stability, (2) availability of another GNSS reference station in the site vicinity, (3) availability of another instrument capable of STD measurements in the site vicinity and (4) the location of the station with respect to its altitude and the weather events which occurred during the evaluation period. The subset also includes collocated (dual) GNSS stations that played an important role in the validation. The collocated stations observed GNSS satellites with the same azimuth and elevation angles, so that they should theoretically deliver the same or very similar tropospheric parameters – ZTD, linear horizontal gradients and slant delays. Post-fit residuals of carrier-phase observations at the collocated stations should represent common effects due to the local tropospheric anisotropy, while systematic differences could remain due to instrumentation and environmental effects such as antenna and receiver characteristics and multipath. Only STDs from the WVR at Potsdam, collocated with the GNSS stations POTM and POTS, were available for this study because the second WVR, located at Lindenberg and collocated with the GNSS stations LDB0 and LDB2, was operated only in the zenith direction during the period of the study. 3 STD retrievals from GNSS observations The STD cannot be estimated directly from GNSS data since the total number of unknown parameters in the solution would be higher than the number of observations. Instead, the total delays in the zenith direction above the GNSS station (i.e. ZTD) are adjusted together with, optionally, total tropospheric linear horizontal gradients (G) to account for the first-order asymmetry of the local troposphere. The estimates are valid for individual processing epochs whenever using a stochastic approach or for a given time interval when modelling the troposphere with a deterministic process, e.g. by a piece-wise constant or linear model. In practice, the ZTD is decomposed into an a priori model, usually by introducing the zenith hydrostatic delay (ZHD; see Saastamoinen, 1972), and the estimated corrections, representing (mainly) the zenith wet delay (ZWD). Similarly, the STD is decomposed to the ZHD, ZWD, Gand post-fit residuals (RES) as described in Eq. (1), where ele is the elevation angle and azi is the azimuth angle in degrees. The STD value is given in metres. STD(ele,azi)=ZHD ·mfh(ele)+ZWD ·mfw(ele) +G(ele,azi)+RES (1) The elevation angle dependency of STD is described by the mapping functions, separately for the hydrostatic (mfh)and the wet (mfw)components. Nowadays, the Vienna Mapping Function (VMF1; see Böhm et al., 2006a) – or VMF1-like concept – is commonly used in GNSS data processing. Also, the empirical mapping function Global Mapping Function (GMF; see Böhm et al., 2006b) is popular since it is consistent with VMF1 and easier to implement (independent on external data needing updates). Both the VMF1 and the GMF are applicable down to 3◦-elevation angles. The first-order horizontally asymmetric delay G(ele, azi) in Eq. (1) reflects local changes in temperature and particularly in water vapour. MacMillan (1995) proposed a model describing the gradient delay as a function of the elevation and azimuth angles: G(ele,azi)=mfg·(GN·cos(azi)+GE·sin(azi)),(2) where mfg(ele)=mfh(ele)·cot(ele). Chen and Herring (1997) replaced the elevation-dependent term mfh(ele)·cot(ele)by the gradient mapping function mfg(ele)=1/(sin(ele)·tan(ele)+C), with C=0.0032, nowadays commonly used in GNSS data processing. Typical range for GNand GEis below 1–2mm, but gradients can reach up to 7mm during extreme weather events. The gradient of 1 mm corresponds to about 55 mm slant delay correction when projected to 7◦-elevation angle. Additionally, post-fit residuals RES may contain unmodelled tropospheric effects not covered by the estimated tropospheric parameters. Such remaining effects are supposed to be caused mainly by higher spatial and temporal variations of the humidity or its significant horizontal asymmetry in the troposphere. Obviously, residuals contain also other un-modelled effects such as multipath, errors in antenna-phase centre variations or satellite clocks. For eliminating such systematic effects, cleaning of post-fit residuals is applied by generating elevationand azimuth-dependent correction maps as described by Shoji et al. (2004). For each solution and each station, we thus computed mean values of post-fit residuals in 1 ×1◦bins using the whole benchmark www.atmos-meas-tech.net/10/2183/2017/ Atmos. Meas. Tech., 10, 2183–2208, 2017 2186 M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays Table 1. Characteristics of 10 GNSS reference stations. Name Latitude Longitude Height Network Dual station Receiver Antenna (◦) (◦) (m) GOPE 49.914 14.786 593 IGS, EPN TPS NET-G3 TPSCR.G3 TPSH KIBG 47.449 12.309 877 TPS GB-1000 TPSCR3_GGD CONE LDB0 52.210 14.118 160 LDB2 JAVAD TRE_G2T JAV_GRANT-G3T NONE LDB2 52.209 14.121 160 LDB0 JPS LEGACY LEIAR25.R4 LEIT POTM 52.379 13.066 145 POTS JAVAD TRE_G3TH JAV_GRANT-G3T NONE POTS 52.379 13.066 144 IGS, EPN POTM JAVAD TRE_G3TH DELTA JAV_RINGANT_G3T NONE SAAL 47.426 12.832 796 TPS GB-1000 TPSCR3_GGD CONE WTZR 49.144 12.879 666 IGS, EPN WTZS, WTZZ LEICA GRX1200+GNSS LEIAR25.R3 LEIT WTZS 49.145 12.895 663 IGS WTZR, WTZZ SEPT POLARX2 LEIAR25.R3 LEIT WTZZ 49.144 12.879 666 IGS WTZR, WTZS JAVAD TRE_G3TH DELTA LEIAR25.R3 LEIT period. Residuals exceeding ±3 times the standard deviation were excluded from the computation of the mean. Computed means were then subtracted from the original post-fit residuals to generate solutions using cleaned residuals. For the analysis of GNSS L1 and L2 carrier-phase observations, the least-squares adjustment or Kalman-filter approach was applied to estimate the ZWDs and the two horizontal gradient components GNand GEat each GNSS site (Table 1). Afterwards, Eq. (1) was used to compute STDs for each satellite in view. Whenever zero-differenced (ZD) post-fit residuals were available for any solution, three variants of the solution are presented in the paper: (1) solution without residuals (nonRES), (2) solution with raw residuals (rawRES) and (3) solution with cleaned residuals (clnRES). Seven institutions delivered their STD solutions for this validation study, namely École Supérieure des Géomètres et Topographes (ESGT CNAM), Geodetic Observatory Pecný (GOP, RIGTC), Helmholtz Centre Potsdam – German Research Centre for Geosciences (GFZ), Royal Observatory of Belgium (ROB), VŠB – Technical University of Ostrava (TUO), Vienna University of Technology (TUW) and Wrocław University of Environmental and Life Sciences (WUELS). Principal information about individual solutions is given in Table 2 with a few specific notes important for the interpretation of the results. GOP delivered two solutions based on the Precise Point Positioning (PPP) technique (Zumberge et al., 1997) and using the in-house developed application Tefnut (Douša and Václavovic, 2014) derived from the G-Nut core library (Václavovic et al., 2013). Considering all available GNSS solutions, only GOP used a stochastic modelling approach to estimate all parameters. Additionally, GOP provided two solutions: (1) GOP_F using Kalman filter (forward filter only), i.e. capable of providing ZTD, tropospheric gradients and STDs in real time; and (2) GOP_S applying the backward smoothing algorithm (Václavovic and Douša, 2015) on top of the Kalman filter in order to improve the quality of all estimated parameters during the batch-processing interval and to avoid effects such as the PPP convergence or re-convergence. Some institutions also delivered two STD solutions which differ in a single processing setting. The aim was to evaluate their impact on STDs: (a) TUO_G and TUO_R exploit GPS-only and GPS +GLONASS observations, respectively; (b) TUW_3 and TUW_7 apply an elevation cut-off angle of 3 and 7◦respectively; and (c) ROB_G and ROB_V use the GMF and VMF1 mapping functions, respectively. Additionally, ROB solutions are the only ones based on the processing of double-difference (DD) observations and providing ZD carrier-phase post-fit residuals converted from the original DD residuals using the technique described in Alber et al. (2000). For other DD solutions, variants without adding residuals were compared only. In total, we validated 11 solutions computed with five different GNSS processing software. Five of the solutions used GPS and GLONASS observations and six solutions used GPS-only observations; five of them are based on DD observations and six of them are computed using zero-difference data in PPP analysis. More information about TUW solutions can be found in Möller et al. (2016), about GFZ in Bender et al. (2009, 2011) and Deng et al. (2011) and about CNAM in Morel et al. (2014). For ROB, TUO and WUE solutions we refer the reader to Dach et al. (2015). 4 Computation of slant total delay from numerical weather prediction model Simulating STDs in NWMs consists in integrating the atmospheric refractivity through the path followed by GNSS signals. STDs have been simulated using three different NWMs: ALADIN-CZ (4.7km resolution limited-area hydrostatic model, operational analysis in 6h interval with forecasts for 0, 1, 2, 3, 5, 6 h; http://www.umr-cnrm.fr/aladin/), ERA-Interim (1◦horizontal resolution, 6h reanalysis) and NCEP-GFS (1◦horizontal resolution, 6 h operational analysis; https://www.ncdc.noaa.gov/data-access/model-data/ model-datasets/global-forcast-system-gfs). None of these NWMs assimilate data from ground GNSS stations. For more information about the models, see Douša et al. (2016) Atmos. Meas. Tech., 10, 2183–2208, 2017 www.atmos-meas-tech.net/10/2183/2017/ M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays 2187 Table 2. Information about individual GNSS-based STD solutions used in the validation. Solution Institution Strategy Software GNSS Elev. Mapping Products ZTD/gradients ZD post-fit name cut-off function interval residuals CNAM ESGT CNAM DD GAMIT GPS 3◦VMF1 IGS final 1 h/1 h No GFZ GFZ Potsdam PPP EPOS 8 GPS 7◦GMF GFZ 15 min/1 h Yes GOP_F GO Pecný PPP G-Nut/Tefnut GPS 7◦GMF IGS final 2.5 min/2.5 min Yes GOP_S GO Pecný PPP G-Nut/Tefnut GPS 7◦GMF IGS final 2.5 min/2.5 min Yes ROB_G ROB DD Bernese 5.2 GPS +GLO 3◦GMF CODE final 15 min/1 h Yes ROB_V ROB DD Bernese 5.2 GPS +GLO 3◦VMF1 CODE final 15 min/1 h Yes TUO_R TU Ostrava DD Bernese 5.2 GPS+GLO 3◦VMF1 CODE final 1 h/3 h No TUO_G TU Ostrava DD Bernese 5.2 GPS 3◦VMF1 CODE final 1 h/3 h No TUW_3 TU Vienna PPP NAPEOS GPS +GLO 3◦GMF ESA final 30 min/1 h Yes TUW_7 TU Vienna PPP NAPEOS GPS +GLO 7◦GMF ESA final 30 min/1 h Yes WUE WUELS PPP Bernese 5.2 GPS 3◦VMF1 CODE final 2.5 min/1 h Yes and specifically Trojáková (2016) for ALADIN-CZ model and Dee et al. (2011) for ERA-Interim. First, STD solutions using the ERA-Interim and NCEP-GFS models were delivered by GFZ Potsdam using the acronym ERA/GFZ and GFS/GFZ, respectively. Only a short introduction is provided in Sect. 4.1 since the GFZ tool for ray tracing has been described in the papers cited below. Two STD solutions were then delivered for the ALADIN-CZ model: (a) the ALA/BIRA, which was generated at the Royal Belgian Institute for Space Aeronomy (BIRA), described in Sect. 4.2; and (b) the ALA/WUELS, which was provided by the Wrocław University of Environmental and Life Sciences, described in Sect. 4.3. 4.1 Description of ERA-Interim STD solution (ERA/GFZ) and NCEP-GPS STD solution (GFS/GFZ) The ERA-Interim and NCEP-GFS STD solutions by GFZ are based on “assembled” STDs. At first, for the considered station and epoch, a set of ray-traced STDs (various elevation and azimuth angles) is computed using technique described in Zus et al. (2014). Secondly, from this set of raytraced STDs, the tropospheric parameters (i.e. zenith delays, mapping function coefficients, firstand higher-order gradient components) are determined. Finally, for the required azimuth and elevation angle the STD is “assembled” using the tropospheric parameters. For a detailed description of the tropospheric parameter determination the reader is referred to Douša et al. (2016). The differences between the “assembled” and ray-traced STDs are sufficiently small in particular for elevation angles above 10◦(Zus et al., 2016). In essence, the largest uncertainty in the “assembled” (and ray-traced) STDs remains the uncertainty of the underlying NWM refractivity field. This uncertainty is estimated to be about 8– 10 mm close to the zenith, increasing to about 8–10cm at an elevation angle of 5◦(Zus et al., 2012). Similar uncertainty of around 8 mm for the zenith direction was also found for ALADIN-CZ model in Douša et al. (2016). 4.2 Description of ALADIN-CZ STD solution from BIRA (ALA/BIRA) To compute STDs from ALADIN-CZ, a simplified strategy has been used to model the curve path followed by GNSS signals through the neutral atmosphere, as suggested by Saastamoinen (1972). The delays simulated with this strategy show small differences in comparison to straight-line simulations (differences of about 4, 5 and 10mm, respectively, at 15, 10 and 5◦elevation). Simulations have computed STDs down to 3◦elevation; however, under an elevation of 15◦a proper ray-tracing strategy, as mentioned in Sect. 4.1, should be applied. For each latitude–longitude grid point and each level of ALADIN-CZ model, the NWM outputs considered to compute STDs are geopotential height (geopotH in m), pressure (Pin Pa), temperature (Tin K), partial pressure of water vapour (ein Pa) and mixing ratio of liquid and solid water (in kg kg−1). The ground pressure of each column is also retrieved. The geopotential height is converted to the altitude above the geoid: hgeoid =(g0·Re·geopotH)/(gc ˙ Re−g0·geopotH), (3) where g0is a standard gravity acceleration (mean value of 9.80665 m s−2from the World Meteorological Organization, WMO); Re=6378137/(1.006803 −0.006706 ×sin2(lat)) is the radius of the ellipsoid in metres for the latitude (lat in degrees); gis the gravity acceleration (in ms−2) of the considered location given as g=9.7803267714 ·(1.+0.00193185138639 ·sin2(lat))/ q1.−0.00669437999013 ·sin2(lat). (4) Then, the height above the geoid is converted to height above the WGS84 ellipsoid (in m) with the use of the EGM96 (Earth Gravitational Model; Lemoine et al., 1998) undulation. Note that for the region of the benchmark campaign the difference between geoid and WGS84 altitude is about 47m. www.atmos-meas-tech.net/10/2183/2017/ Atmos. Meas. Tech., 10, 2183–2208, 2017 2188 M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays Using the hypsometric equation, the ground pressure and the pressure of each level are considered to estimate the altitude for the different levels. In total ALADIN-CZ outputs provide 87 levels up to an altitude of about 55km. However, to assess STDs from ALADIN-CZ, the integration was stopped at 15km since the contribution of water vapour above this altitude is negligible. An adaptive step is considered (100, 200, 250, 500, 1000 m, respectively, for vertical altitudes between 0 and 1, 1 and 3, 3 and 5, 5 and 10, and 10 and 15km). Bi-linear interpolations of ALADIN-CZ parameters at the altitude of the GNSS station and for each step of the integration were proceeded. Note that there is no station selected for the validation located below the first layer of ALADIN-CZ. The expression of simulated STDs from ALADIN-CZ is the summation of these four contributions: STD =SHDint +SWDint +SHMDint +STDext,(5) where SHDint, SWDint and SHMDint are, respectively, the inside-model integration contribution of the hydrostatic, wet and hydrometeor delays, and STDext is the external model contribution (over 15km). SHDint =10−6Xk=ktop k=1k1Pi Tvi 1si,(6) SWDint =10−6Xk=ktop k=1 k02ei Ti +k3ei T2 i!1si,(7) with k0 2=k2−k1·Rd/Rw, where k1(in K Pa−1), k2(in K Pa−1) and k3(in K2Pa−1) are the empirical refractivity coefficients of Bevis et al. (1994); Rdand Rwthe gas constants, respectively, for dry air and water vapour (in J kmol−1K−1); and Tvis the virtual temperature (in K). For the estimation of the hydrometeor contribution inside the model, as presented in Eq. (8), (Nlw,Mlw)and (Nice,Mice)are atmospheric refractivity and mass content per unit of air volume of the liquid and ice water, respectively. SHMDint =10−6Xk=ktop k=1(Nlw +Nice)1si =Xk=ktop k=1(αlwMlw +αiceMice)1si(8) The estimation of coefficients αlw ∼1.45 and αice ∼0.69 is presented in Brenot (2006). The ALADIN-CZ model provides mixing ratios of cloud water (liquid components) and pristine ice (solid water components). The mass content per unit of air volume is obtained using the associated mixing ratio, pressure, water vapour partial pressure and temperature. STDext is obtained with the hydrostatic formulation (Saastamoinen, 1972) mapped with mfh(see Eq. 1) and using the elevation, latitude and pressure of the last step of the integration (i.e. at 15 km). Note that the wet contribution over 15 km is neglected since it is practically zero. The estimation of STDext (about 0.21m) provides sufficiently accurate modelling for the hydrostatic contribution over 15 km (as shown by the sensitivity test from Brenot, 2006). 4.3 Description of ALADIN-CZ solution from WUELS (ALA/WUELS) The ray-traced tropospheric delays for WUELS’ solution are based on piece-wise bent-2-D model propagation. Thus, it prevents us from knowing the exact trajectory in advance, in contrast to straight-line models, and must be solved iteratively based on the preceding ray refractive index. Similar examples are given by Böhm and Schuh (2003) and Hobiger et al. (2008). We assume the ray path does not leave the plane of constant azimuth for a given elevation angle to a satellite. The out-of-plane contribution to the delay is thus neglected, making the propagation two-dimensional (hence 2-D). The real ray path is then approximated by a finite number of linear ray pieces in WGS84 coordinates using Euler’s formula for the Earth radius: R=(cos2A/M +sin2A/N)−1,(9) where Ais the azimuth angle between a satellite and a receiver, and Mand Nare radii of curvature along meridian and prime vertical, respectively. We follow height-dependent increments as presented in Rocken et al. (2001): 10, 20, 50, 100 and 500 m, respectively, for geometric altitudes between 0 and 2, 2 and 6, 6 and 16, and 16 and 36km and above 36 km, which require meteorological parameters to be vertically interpolated in order to obtain finer resolution. Both Pand eare interpolated exponentially from the two nearest layers, while the temperature change is considered linear. Horizontally, we find the four nearest nodes for each ray to perform weighted mean interpolation, where the weighting function equals the inverse squared distance. The reference hybrid level of the ALADIN-CZ model is determined by surface geopotential, which is converted to geopotential metres by dividing the geopotential values by g0. Meteorological parameters are expressed on pressure levels which represent standard vertical coordinates. The hypsometric equation is used to calculate geometric thickness between consecutive isobaric surfaces: dz=Rd·Tm/g0·ln(P1/P2), (10) where Rd=287.058 JK−1kg−1is the gas constant for dry air, and Tmis the mean virtual temperature of the layer between P1and P2pressure levels in Kelvin. The conversion from ALADIN-CZ vertical coordinate system to geometric altitudes is then consistent with the BIRA approach described in Sect. 4.2. In WUELS’ solution, the signal tracking is performed exploiting a full model vertical resolution up to the uppermost ALADIN-CZ layer at 55km. Above the top layer, we adopt the US Standard Atmosphere (1976) to provide supplementary meteorological data up to 86km. For each ray-path coordinate, the refractive index is calculated as a function of P(in hPa), e(in hPa) and T(in K) with empirically derived “best available” coefficients kgiven by Rueger (2002). N=(n−1)×106=k1·(P −e)/T +k2·e/T +k3·e/T 2(11) Atmos. Meas. Tech., 10, 2183–2208, 2017 www.atmos-meas-tech.net/10/2183/2017/ M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays 2189 The contribution of water droplets and ice crystals in the atmosphere is neglected in the total delay. All tropospheric delays are traced with respect to vacuum elevation angles. The electromagnetic delay is calculated for a given chord length (s) using the mean refractive index nbetween two consecutive rays, yielding the total delay in metres: STD =Xsi(ni−1)×106,(12) which can be separated on hydrostatic and wet part using respective refractive indices. Additionally to the radio path length, the accumulated bending effect (bend) along the ray path is added to the hydrostatic mapping function which, together with the wet mapping function, can be calculated as follows: bend =X(si−cos(elei−elek)si), (13) mfh=(SHD +bend)/ZHD,(14) mfw=SWD/ZWD,(15) where eleiis the elevation angle for a given model layer and elekis the outgoing elevation angle at uppermost altitude. 4.4 Assessment of the hydrostatic, wet and hydrometeor contributions to the slant delays ALADIN-CZ NWM has been used to estimate the hydrostatic, wet and hydrometeor contributions to slant delays. During the whole period of the benchmark campaign, the maximum contribution of hydrometeors reached 17mm at the zenith during the extreme weather events on 20–23 June (Douša et al., 2016). The 2-D fields of ZTD, ZHD, ZWD and ZHMD (zenith hydrometeor delays) are presented in Fig. 1. They illustrate the large-scale convection with the presence of hydrometeors along the convergence line associated with a strong contrast of dry and wet air masses. The contribution of hydrometeors to ZTD reached up to 7mm (as scaled in the zenith direction) for the stations POTS and POTM at 15:00 UTC on 23 June 2013 (see Fig. 1d). According to satellite trajectories at this time for the station POTS, a maximum SHMD of 25.6mm is observed for a satellite at 22◦-elevation angle. Figure 2 shows simulated differential STDs for a cone with a 10◦-elevation angle during the severe weather condition of the 23 June 2013 and mapped in the zenith direction (at 90◦) using the mapping functions of Eq. (1): mfh for SHD and mfwfor SWD and SHMD. For this 10◦cone, the minimum present values of total, hydrostatic, wet and hydrometeors delays simulated at 15:00UTC on 23 June are given as STDmin, SHDmin, SWDmin and SHMDmin in Fig. 2. The respective differences of STD, SHD, SWD and SHMD and corresponding minimum values simulated at 15:00 UTC are presented in Fig. 2. The anisotropic variation of total, hydrostatic, wet and hydrometeor delays can be visualized on a skyplot. As a confirmation of Figs. 1b and 1d, 2 shows weak hydrostatic anisotropy. This anisotropy (up to 5.8 mm) is almost the same as the hydrometeors one (up to 6 mm). The area within the red curve is larger than the purple area (hydrostatic anisotropy), meaning that the total effect of the hydrometeor anisotropy is slightly larger than the one from the hydrostatic component. Note that Fig. 2 shows the anisotropies simulated at 10◦and mapped at 90◦ (giving an idea of the variations in the zenith direction). The largest anisotropy is clearly induced by water vapour (values up to 20 mm in the south-east direction of POTS, also shown in Fig. 1c). With mean hydrostatic and hydrometeor anisotropies oriented in the opposite direction of the wet one, Fig. 2 presents a total anisotropy with weaker values (up to 12 mm) than the wet anisotropy. To complement the snapshot of Fig. 2 the time evolutions of SHD, SWD, SHMD and STD in the direction of all observed GNSS satellites for the station POTS are presented in Fig. 3. Slant delays have been simulated in the direction of observed satellites (hydrostatic, wet, hydrometeor and total contributions) and, to avoid the effect of the elevation and to look at the same order of magnitude of delays, corresponding delays in the zenith direction have been computed and mapped using mapping functions presented in Eq. (1) (mfhfor SHD and mfwfor SWD and SHMD). These values (STDmin =2310.6 mm, SHDmin =2240.8 mm, SWDmin =43.1 mm and SHMDmin =0 mm), obtained during the whole period of the benchmark campaign, have been subtracted from their corresponding values simulated in direction of satellites. Then, the differences have been mapped back at 90◦. For day of year (DOY) 174 (i.e. 23 June 2013), we can see a contribution of hydrometeors up to 10mm. Looking at the whole period of the benchmark campaign, the variation ranges of STD, SHD and SWD (mapped at 90◦) are 275, 80 and 230 mm, respectively. Figure 3 illustrates the use of GNSS delay observations in meteorology (detection of variation of water vapour represented by the wet delay, as well as detection of heterogeneities from hydrostatic delays and occasionally from hydrometeors in specific severe weather cases). Note that there are no data available for POTS station between DOY 121 and 125 and DOY 160 and 163. For this reason, we have not simulated the slant delays for this period, as shown by the gaps in Fig. 3. The simplified strategy used to simulate curve slant paths gives some inaccurate simulations of slant delays for elevations between 3 and 5◦, shown by isolated values in Fig. 3. Such inaccuracies could be avoided by using a ray-tracing algorithm. For a comprehensive overview on ray-tracing algorithms and comparisons the reader is referred to Nafisi et al. (2012). 5 Water vapour radiometer measurements During the benchmark period, the WVR located at GFZ Potsdam operated in a mode that scanned the atmosphere at selected elevation and azimuth angles. The instrument is situated on the same roof as the GNSS reference stations POTM www.atmos-meas-tech.net/10/2183/2017/ Atmos. Meas. Tech., 10, 2183–2208, 2017 2190 M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays Figure 1. Simulation of ZTD, ZHD, ZWD and ZHMD at 15:00UTC on 23 June 2013. Each black dot represents a GNSS station included in the benchmark dataset. For stations included in this STD validation study their names are given. and POTS. All three devices are within 10m from each other. The HATPRO WVR from Radiometer Physics was set up to scan the atmosphere to extract profiles of atmospheric temperature, water vapour and liquid water using frequencies between 22.24 and 27.84GHz and a window channel at 31.4 GHz. The WVR switches between “zenith mode” when it measures IWV and “slant mode” when it tracks GPS satellites using an in-built GPS receiver. In the latter case, SIWV values are delivered for the direction of satellites. Since the instrument can track only one satellite at one moment the number of observations is quite limited compared to slants from GNSS that are simultaneously observed from several GNSS satellites. Our study focuses on the comparison of STDs, not SIWV. It was thus necessary to convert the WVR SIWV into STDs. Firstly, WVR observations with rain flag and atmospheric liquid water (ALW) values exceeding 1kg m−2were rejected. Both rain and high values of ALW can significantly distort the quality of WVR measurements. Secondly, SIWV values were converted into SWDs using the Askne and Nordius (1987) formula and the refractivity constants from Bevis et al. (1994). ZHD values were computed with the precise model given by Saastamoinen (1972). For the described conversions, we used values of the atmospheric pressure and temperature measured in situ of the GNSS reference station POTS. A hydrostatic correction for the altitude difference between the meteorological station and the WVR position was applied to the atmospheric pressure values. ZHD values were mapped to elevation angles of the WVR using the hydrostatic mapping function derived from the NCEP-GFS (Douša et al., Atmos. Meas. Tech., 10, 2183–2208, 2017 www.atmos-meas-tech.net/10/2183/2017/ M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays 2191 Figure 2. Skyplot of differential slant delays simulated at 10◦ and mapped at 90◦, for a 360◦azimuthal range (at 15:00UTC on 23 June 2013). For total, hydrostatic, wet and hydrometeors delays, a differential slant delay is the difference between a slant delay simulated and the respective minimum value (obtained considering slant delays simulated at 10◦elevation along all the azimuthal directions). Figure 3. Time series of slant delays (STD, SHD, SWD and SHMD) differences (in direction of all GNSS visible satellites, then mapped in the zenith direction) during the whole period of the benchmark campaign for the station POTS. 2016). In order to convert accurately SIWV to STDs, we took into account the influence of the hydrostatic horizontal gradients (see e.g. Li et al., 2015b). We used the hydrostatic horizontal gradients derived from the NCEP-GFS for that purpose. Finally, SHD and SWD values were summed up to deliver STDs. The described conversion of WVR SIWV to STDs aimed to minimally distort the accuracy of original WVR observations. 6 Methodology of STD comparisons We provide the specificities of each type of technique comparisons in this section. Since NWM outputs are restricted to the time resolution of their predictions (typically 1, 3 or 6 h) and, since WVR is able to track only one satellite at one moment, all three sources provide different numbers of STDs per day. Therefore, three different comparisons are presented: (1) results for GNSS versus GNSS comparisons, (2) results for GNSS versus NWM comparisons and (3) results for GNSS versus WVR comparisons. Section 7 presents the validation at individual stations and Sect. 8 intercompares results obtained at GNSS dual stations. All the given results are obtained over the whole benchmark period. No outlier detection and removal procedure was applied during the statistics computation within the study. Two variants of the comparisons are presented: “ZENITH” and “SLANT”. “ZENITH” stands for original STDs mapped back to zenith direction using 1/sin(e) formula. Such mapping aimed to normalize STD differences for their evaluation in a single unit. The “SLANT” type of comparison denotes an evaluation of STDs at their actual elevation angles. To be more specific, slant delays were grouped into individual elevation bins of 5◦; i.e. for example all slants with an elevation angle between 10 and 15◦were evaluated as a single unit. There was one exception regarding the size of a bin since the lowest one contained slants from elevation angles of 7 to 10◦, 7◦being the lowest elevation angle common to all GNSS STD solutions. This cut-off angle was thus used in all GNSS versus GNSS and GNSS versus NWM comparisons. Presented values of biases and standard deviations were computed directly from all STDs within the processed benchmark campaign period, and therefore they are not based on any kind of daily or other averaging. In some tables, only median values of bias and standard deviation over all GNSS STD solutions (Tables 5, 7 and 8) or over all processed stations (Tables 3 and 4) are given to consolidate the presentation of validation results. Median was used as a parameter minimally affected by outliers. 6.1 GNSS versus GNSS comparisons In the case of individual inter-GNSS solutions validation, the situation was straightforward and no interpolation nor specific hypothesis was necessary: the comparisons were done on a direct point-to-point basis of observations coming from identical azimuth and elevation directions. www.atmos-meas-tech.net/10/2183/2017/ Atmos. Meas. Tech., 10, 2183–2208, 2017 2198 M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays Figure 8. Comparison of individual GNSS STD solutions without residuals (nonRES) against NWM solutions ALA/BIRA, ERA/GFZ, GFS/GFZ and ALA/WUELS (from top to bottom), projected in the zenith direction: bias (a, c, e, g) and standard deviation (b, d, f, h). good agreement between these independent sources used for retrieving slant delays. Standard deviations generally range from 8 to 12mm in the zenith projection, with the exception of ALA/WUELS, which shows lower precision by a factor of 2.5. Statistics stem from the complete benchmark period, and it should be noted that the daily variation of GNSS STDs was much lower than of NWM ray-traced STDs. Significantly higher values of biases and standard deviations were observed at particular days for NWM solutions. A detail evaluation of daily statistics with respect to the extreme weather conditions is one of the topics that we will study in future. 7.3 GNSS versus WVR Figure 10 compares GNSS and WVR solutions at stations POTM and POTS, in the zenith direction. The number of Atmos. Meas. Tech., 10, 2183–2208, 2017 www.atmos-meas-tech.net/10/2183/2017/ M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays 2199 Figure 9. Comparison of NWM-based solutions (ALA/BIRA, ERA/GFZ and GFS/GFZ) against GNSS GFZ solution at station POTS, in the slant direction. slant observations which entered the comparison was 32794 at station POTM and 36 070 at station POTS. Two remarks can be made on the evaluation of biases. Firstly, an overall bias of about 4mm between the stations POTM and POTS, visible for all GNSS solutions already in Fig. 8, indicates a common issue with the GNSS data processing at the station POTM. It was particularly increased for GOP_F, GOP_S and GFZ PPP solutions. Secondly, a bias of about 5.5 mm in the zenith direction can be found between WVR and GNSS solutions even at station POTS. This bias roughly corresponds to 1 kg m−2of IWV, which can be considered the achievable accuracy of either of the two techniques; however, WVR accuracy is more dependent on a proper instrument calibration. Values of standard deviation, resulting mostly in 12 mm, are higher than those observed in GNSS versus GNSS comparisons (Sect. 7.1) and slightly higher than from GNSS versus NWM comparisons (Sect. 7.2). A cut-off elevation angle of 15◦was used for the comparison with WVR STDs instead of 7◦used in other validations. Additionally, it has to be noted that the results can be partly influenced with the settings applied for finding pairs between GNSS and WVR STDs (Sect. 6). STDs from WVR can thus originate from slightly different azimuth/elevation angles and times than the GNSS ones. All GNSS solutions perform similarly against WVR, with the exception of GOP_F due to the application of a real-time capable strategy. The GNSS versus WVR validation at the station POTS using original elevation angles is displayed in Fig. 11. Although Figure 10. Comparison of individual GNSS STD solutions for stations POTM and POTS versus WVR measurements, expressed in the zenith direction, bias (a) and standard deviation (b). The median value of all solutions at each station is represented by the dotted blue line in each bin. some differences between GNSS solutions are visible, all of them performed in a very similar manner. The decrease of values of four statistical parameters strongly follows the increase of elevation angle and, generally, it is steeper than statistics dependency of GNSS versus NWM. It indicates that slant delays from WVR below 40◦become generally unreliable, which is particularly clear from normalized biases and www.atmos-meas-tech.net/10/2183/2017/ Atmos. Meas. Tech., 10, 2183–2208, 2017 2200 M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays standard deviations at lower-elevation angles. A sudden increase of the values is observed at elevations of 55–60◦, most likely originating from WVR observations which are not yet understood. Generally, standard deviations for all solutions using cleaned residuals (raw residuals) are on average 1.7% (3.8 %) higher than for the solutions without residuals. Differences between solutions variants are smaller due to an overall higher uncertainty of WVR observations, but the results are in a good agreement with those obtained for GNSS versus NWM comparisons presented in Sect. 7.1. 8 Validation of results at collocated stations Two erroneous techniques for STD retrievals have been compared in previous sections (GNSS vs. NWM, GNSS vs. WVR) without knowing the true reference. The errors stem from the observation noise on one hand and from the processing models including the model for adjusted parameters on the other hand. From this perspective, the higher standard deviations for GNSS STD solutions applying clean residuals compared to those using adjusted GNSS parameters only (without residuals) do not necessarily mean the lower quality of the former. GNSS and NWM models with limited temporal and spatial approximations are not able to represent true signal tropospheric delays between a receiver and all visible satellites. The simplifications certainly result in better agreement of STDs without residuals in Eq. (1), but they hardly represent the true tropospheric path delays, deviating particularly during the events with high spatiotemporal variations in the troposphere. For this reason, we assessed all GNSS solutions at the collocated (dual) stations because for such constellation we are able to provide troposphere-free differences of STDs to evaluate noise of GNSS STD retrievals. We particularly focused on days with a high variability in the troposphere selected from the benchmark period. Dual stations were available in the benchmark campaign at three different locations in Germany. The first two sites collocate twin GNSS reference stations (LDB0 +LDB2 and POTM +POTS) and the third location collocates three individual reference stations (WTZR +WTZS +WTZZ). Nevertheless, in the case of Wettzell, only results for WTZR+WTZS are presented due to their similarity with the two other combinations at the same place. Characteristics of the stations are summarized in Table 6. 8.1 Slant residuals and slant delay differences STD validations in this paper were done for 2 months of the benchmark period during which heavy rain events occurred for some days, particularly 31 May–3 June, 9–11 and 21– 26 June, all causing severe flooding in central Europe. During normal weather conditions, the tropospheric variation is reasonably smooth, meaning it can be well represented by GNSS STDs reconstructed from ZTDs and horizontal gradients. However, during high temporal or spatial variabilities in the troposphere, post-fit residuals certainly contain tropospheric signals which were not modelled. If they surpass the observation noise and other residual errors from GNSS models, cleaned residuals should be considered in the GNSS STD model as described in Eq. (1). In order to initially address optimal STD modelling under different weather conditions within the benchmark period, we tried to identify days with a high variability in the troposphere. Daily standard deviations of cleaned post-fit residuals were computed individually for each day of the benchmark period, for every station and GNSS solution for 1◦-elevationangle bins. We studied their daily variations considering the GNSS model applied. If cleaned post-fit residuals consist of the noise of observations only, the variation in time should be negligible. However, the days showing significantly higher values, correlated at all collocated stations, indicate highly variable tropospheric conditions. Three such days were identified at LDB0, LDB2, POTM and POTS stations (31 May, 20 June, 23 June) and 2 days at WTZR and WTZS stations (19 and 20 June). They all very well correspond to the days initiating heavy precipitations in the domain (Douša et al., 2016). Typical differences between raw and clean residuals are displayed in Fig. 12 for all elevations during the normal day (19 June, DOY 170) and the day with high variability in the troposphere (DOY 171, 20 June) for LDB0, LDB2, POTM and POTS stations using GFZ solution. Obviously, the variability of clean residuals (black dots) and their 2σenvelopes are higher by a factor of 2 for the day of year 171 compared to 170. The variability is clearly visible over all elevations, but the increase is slightly higher at low elevations. The plots for these four stations clearly demonstrate the different quality of GNSS observations, particularly related to a multipath effect displayed by 2σenvelope (green curves). A low multipath is common to the stations using choke ring antennas, in our case POTS and LDB2, but LDB2 still suffers from unknown systematic effects at 35–55◦elevations. A very high multipath effect was observed at LDB0 station over all elevations. Variability of 2σenvelopes of clean residuals (red curves) indicates a higher sensitivity of clean residuals to the weather conditions compared to station selection and observation quality, thus suggesting a significant contribution from the troposphere to the cleaned residuals. In the same context, raw residuals show much higher sensitivity to the observation quality compared to different weather conditions, which is particularly true in the case of LDB0 and LDB2 stations. Elevation-dependent differences of STDs using clean residuals (black dots) are displayed in Fig. 13 for the same days as in Fig. 12, selecting GFZ solution and station pairs WTZS–WTZR and LDB0–LDB2. Additionally, 2σenvelopes are plotted for differences without residuals (red Atmos. Meas. Tech., 10, 2183–2208, 2017 www.atmos-meas-tech.net/10/2183/2017/ M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays 2201 Figure 11. Comparison of WVR against individual GNSS STD solutions at station POTS, in the slant direction. Table 6. Characteristics of individual dual stations. Dual station Location Horizontal Vertical Identical Identical Pairs of distance distance type of type of observations (m) (m) receiver antenna LDB0 +LDB2 Lindenberg 177 0.6 No No 143 005 POTM +POTS Potsdam 2.5 −0.5 No No 180 636 WTZR +WTZS Wettzell 69 2.6 No Yes 84 443 curves), clean residuals (green curves) and raw residuals (blue curves). Firstly, we note that STD differences are more or less similar for both days, i.e. not significantly different between days with normal and high variations in the troposphere, which is also found for other days of the benchmark period. It suggests that increased residuals in Fig. 12 for DOY 171 contain strong contributions from the tropospheric effect that could not have been assimilated into ZTDs and tropospheric horizontal gradients. An alternative explanation suggests a possible contribution of satellite-specific errors common to both receivers, thus easily eliminated in STD differences at the dual stations. However, systematic errors at satellites are well absorbed by initial phase ambiguities in PPP and short-term or random errors, e.g. due to satellite clocks, are in this study eliminated by the use of final products, i.e. stable enough to avoid observed day-to-day variability in cleaned residuals. The DOY 171 thus shows the situation when cleaned residuals contain a tropospheric signal that should be added to the STD retrievals. In the case of GFZ, the contribution from residuals is particularly important due to local troposphere variation in time when using model of piece-wise constant function with 15min time resolution for ZTD and 60 min for horizontal gradients. It is not so obvious in the case of a stochastic process used for epoch-wise estimates of all tropospheric parameters. However, the uncertainty of estimated parameters is then higher compared to the deterministic model, which makes it more difficult to separate errors in estimated parameter and errors due to insufficiency of the linearized tropospheric model in time. Secondly, we can see that envelopes of differences using raw residuals are always the largest ones. Raw residuals vary more with the elevation angle, which is particularly visible for differences between LDB0 and LDB2. Obviously, it is due to the large systematic errors at LDB0 station and additional contribution from LDB2 errors observed at 35– www.atmos-meas-tech.net/10/2183/2017/ Atmos. Meas. Tech., 10, 2183–2208, 2017 2202 M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays Table 7. Comparison of GNSS STDs from the elevation angles ranging from 7 to 15◦at three dual stations; results for days with high daily variability of cleaned post-fit residuals (top) and results for days with low daily variability of post-fit residuals (bottom). Median values of biases and standard deviations (SD) calculated over all GNSS STD solutions are given; statistics are expressed in the zenith direction. nonRES clnRES rawRES Bias (mm) SD (mm) Bias (mm) SD (mm) Bias (mm) SD (mm) Days with high variability of post-fit residuals LDB0 +LDB2 −1.56 4.63 −1.44 5.51 −1.52 5.89 POTM +POTS −5.24 1.89 −5.16 3.47 −5.91 4.24 WTZR +WTZS −0.24 2.31 −0.06 3.25 −0.03 3.77 Days with low variability of post-fit residuals LDB0 +LDB2 −0.52 3.06 −0.52 4.23 −0.59 5.05 POTM +POTS −4.97 1.87 −5.03 3.00 −5.79 3.87 WTZR +WTZS −0.01 1.87 −0.05 3.22 −0.09 3.82 Table 8. Comparison of GNSS STDs from the elevation angles ranging from 15 to 90◦at three dual stations; results for days with high daily variability of cleaned post-fit residuals (top) and results for days with low daily variability of post-fit residuals (bottom). Median values of biases and standard deviations (SD) calculated over all GNSS STD solutions are given; statistics are expressed in the zenith direction. nonRES clnRES rawRES Bias (mm) SD (mm) Bias (mm) SD (mm) Bias (mm) SD (mm) Days with high variability of post-fit residuals LDB0 +LDB2 −0.92 4.19 −0.88 6.13 −0.88 8.87 POTM +POTS −5.28 1.68 −5.30 3.39 −5.30 4.64 WTZR +WTZS −0.53 2.29 −0.56 4.12 −0.49 5.21 Days with low variability of post-fit residuals LDB0 +LDB2 −0.07 2.59 −0.07 4.93 0.04 7.83 POTM +POTS −4.96 1.82 −4.92 3.30 −4.93 4.65 WTZR +WTZS −0.01 1.74 −0.05 3.77 0.02 5.07 55◦elevations. The 2σenvelopes of STD differences with clean residuals smoothly follows the 2σenvelope of STDs differences without residuals, keeping the difference within ±15 mm over all elevations. This indicates a stable and reliable usage of clean residuals under any conditions. In contrast, applying raw residuals at problematic sites may seriously degrade STDs as observed at LDB0 station. Finally, we can consider error contribution from both stations to STD differences at dual stations equal, i.e. δ2 STD_dif =2δ2 STD_res,(18) with δ2 STD_dif variance calculated from cleaned STD differences at specific elevations when using the same processing strategy at both dual stations and with δ2 STD_res characterizing the variance over errors in GNSS STD retrievals corresponding to the observation elevation angle and the applied strategy. Although we can note some differences in δ2 STD_dif in collocations, partly due to differences in contributions from both stations, the relative performance of differences from STDs with clean residuals (green curves) and without residuals (red curves) for different days remains similar. Uncertainties of the simplified STDs at low elevations surpass additional uncertainties due to applying clean residuals (green curves vs. red curves). According to the magnitude of clean residuals at these elevations (Fig. 12), the small uncertainties from calculated differences indicate the presence of tropospheric signals in the residuals at low elevations, roughly below 30◦. It seems to be almost independent from the weather conditions and is supposed to represent mainly unmodelled horizontal asymmetry in the troposphere. However, further study on detail impact of residuals on GNSS STDs modelling during severe weather conditions requires longer datasets, which will be subject of our upcoming study. Figure 14 displays results for comparisons of individual dual stations in slant directions calculated from all days of the benchmark period. The same statistics and plots (not displayed) were prepared also for days identified with “severe” weather conditions, but only minor differences were observed. Strong variations are observed mainly in normalAtmos. Meas. Tech., 10, 2183–2208, 2017 www.atmos-meas-tech.net/10/2183/2017/ M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays 2203 Figure 12. Elevation-dependent variability of clean residuals (black dots) and their 2σenvelopes (red curves) are showed for 19 June (DOY 171) and 20 June (DOY 170) and four stations: POTS, POTM, LDB0 and LDB2. Additionally, plots display 2σenvelopes for raw residuals (blue curves) and multipath (green curves). ized biases over all elevation angles for the solutions using raw post-fit residuals (rawRES) regardless weather conditions. These are clearly related to local effects such as multipath or modelling instrument-related effects (phase centre offsets and variations) and disappear after using the cleaned residuals (clnRES). The standard deviations and normalized standard deviations at all stations are clearly the lowest for variants without using post-fit residuals (nonRES), slightly higher using cleaned residuals and significantly higher when using raw residuals, i.e. corresponding to above-performed inter-technique validations. 8.2 Differences in zenith direction Tables 7 and 8 show the statistics expressed in the zenith direction for observations ranging in elevation angles from 7 to 15 and from 15 to 90◦, respectively. Median values computed over all GNSS solutions for which residuals were available are presented. Results for the identified days with high daily tropospheric variation are given in the upper part of the table – days are stated in the previous section. In the bottom part results for selected days with low daily variation of post-fit residuals are presented. These days were the same for all collocated stations: 25 May, 30 May and 6 June (DOY 145, 150, 157). Biases remain stable regardless of the severe weather occurrence and whether post-fit residuals are used. The lowest standard deviations for all dual stations are always related to the solutions without using post-fit residuals. Interestingly, when comparing the statistics for STDs evaluated separately for ranges of 7–15 and 15–90 elevation degrees, standard deviations are smaller in high compared to low elevations for variants without using residuals and vice versa for variants using either cleaned or raw residuals. This can be interpreted as the standard GNSS tropospheric model (ZTD and horizontal gradients) representing well observations at elevawww.atmos-meas-tech.net/10/2183/2017/ Atmos. Meas. Tech., 10, 2183–2208, 2017 2204 M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays Figure 13. Elevation-dependent variability in STD differences of clean residuals (black dots) and their 2σenvelopes (green curves) are showed for 19 June (DOY 171) and 20 June (DOY 170) and two dual-stations: WTZS–WTZR and LDB0–LDB2. Plots also display 2σ envelopes for differences of raw residuals (blue curves) and without residuals (red curves). tions above 15◦but suffering under the modelling deficiencies mainly at low elevations. These statistics also support the above statement that cleaned residuals are valuable particularly for reconstructing low-elevation STDs regardless of the weather conditions as they certainly contain non-negligible tropospheric signals from high-order horizontal asymmetry. During the days with high variation of residuals, the standard deviations are usually a little bit higher than during the days with low variation, but there is no difference between these two regarding the above-mentioned behaviour. 9 Conclusions We presented results of validating tropospheric slant total delays obtained from GNSS data processing with those obtained from NWM ray tracing, WVR measurements and collocated GNSS stations, in search of the optimal method for estimating GNSS STDs. Ten GNSS reference stations were selected, exploiting data from the 56-day COST ES1206 benchmark campaign. Eleven GNSS solutions, four NWMbased solutions and one WVR-based dataset entered this validation study. Eight out of 11 GNSS solutions delivered STDs in three variants: (1) without post-fit residuals, (2) with raw post-fit residuals and (3) with cleaned post-fit residuals. The comparisons were carried out into two scenarios, firstly for STDs at their true elevation angles and secondly for STD differences mapped into the zenith direction using a simple mapping function. Comparisons of STD solutions without residuals, with raw or with cleaned residuals were used to study the impact of different strategies for optimally retrieving STDs from GNSS. The impact of cleaning residuals led to the standard deviations reduced by a factor of 1.2–1.5 over all stations and solutions, namely reaching 2.5–4.5mm in the zenith direction for clean residuals compared to 3.0–6.5 mm for raw residuals, the latter also being highly dependent on the station. The impact of adding raw or cleaned residuals was practically negligible in terms of biases, which always remained within ±0.1 mm. Biases and standard deviations between GNSS and NWM solutions depended on applied ray-tracing method, NWM source and station location. Worse results, by a factor of 2.5 in terms of standard deviation, were observed for the ALA/WUELS solution originating from the deficiency of the applied ray-tracing method. Generally, biases in the zenith direction were below ±3mm for other solutions with the exception of a positive bias of 5mm observed for GFS NWM model. Standard deviations for all GNSS versus NWM STD comparisons were at the level of 10 mm, excluding the ALA/WUELS solution. Contrary to the GNSS versus GNSS comparisons, normalized standard deviations showed pronounced variability with the elevation angle. Using the simulation of delays from ALADIN-CZ weather model, we illustrated the impact of the hydrostatic, wet and hydrometeors contributions to zenith and slant delays. These showed strong horizontal variations that allowed relevant characterization of mesoscale meteorological situations. Visualizing the slant anisotropic variation of total, hydrostatic, wet and hydrometeor delays in a common skyplot illustrated a weak hydrostatic anisotropy (up to 5.8mm) that was almost the same as that of the hydrometeor (up to 6mm). The largest anisotropy was induced by water vapour (up to 20 mm), but the total anisotropy was much weaker (12mm) due to the compensation of mean hydrostatic and hydrometeor anisotropies oriented in the opposite direction. Atmos. Meas. Tech., 10, 2183–2208, 2017 www.atmos-meas-tech.net/10/2183/2017/ M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays 2205 Figure 14. Comparison of GNSS STDs at dual stations computed over whole benchmark period from individual GNSS solutions in the slant direction for dual stations from left to right: LDB0-LDB2, POTM-POTS and WTZR-WTZS. Statistical parameters from top to bottom: bias, normalized bias, standard deviation and normalized standard deviation. GNSS STDs from stations POTM and POTS were validated against collocated WVR observations pointed to GNSS satellites. A positive bias of about 5.5 and 10 mm was observed for POTS and POTM station, respectively. Standard deviations from comparisons of GNSS versus WVR STDs reached 12 mm in the zenith direction, thus higher compare to NWM solutions. Normalized standard deviations revealed a strong elevation dependency, indicating the WVR observations lack this quality at low elevations, particularly below 40◦. Collocated GNSS stations at three different locations were used to evaluate the quality of GNSS STD retrievals applying statistics over troposphere-free STD differences from theoretical point of view. We could observe strong systematic errors in raw residuals at any elevation angles, particularly at stations without the choke ring antenna, such as LDB0 and POTM. We found a strong elevation dependency of bias when using raw residuals which almost vanished when cleaning the residuals from visible systematic errors. This suggests that the use of raw residuals should not be recommended, at least not without any information about possible systematic errors. Although the simplified STDs reconstructed from the estimated GNSS tropospheric parameters performed the best in all the comparisons, it obviously missed part of tropospheric signals due to non-linear temporal and spatial variations in the troposphere. By identifying low and high variability in the troposphere during all days in the benchmark period, we showed that residuals contain significant tropospheric signals in addition to the simplified model, particularly during high variability in the troposphere. Additionally, we also identified tropospheric signals at low elevations due to a non-linear horizontal asymmetry in cleaned residuals regardless of the station selection and the quality of its observations. From such findings, we recommend the use of cleaned residuals for an optimal STD retrievals from GNSS, at least for low-elevation angles and during high variability in the troposphere. We also have not seen any obvious degradation of STD retrievals in other conditions. The better inter-solution and inter-technique agreements of STDs without residuals compared to those using clean www.atmos-meas-tech.net/10/2183/2017/ Atmos. Meas. Tech., 10, 2183–2208, 2017 2206 M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays residuals are attributed to the too-simple tropospheric model resulting in smooth and robust STDs and, consequently, not containing all interesting signals from the troposphere. The majority of evaluated GNSS solutions used deterministic models with rather long validity of estimated tropospheric parameters for which the residuals are important to overcome modelling deficiencies of low-resolution parameter estimates in time. Our future study will focus on the evaluation of GNSS STDs estimated using a stochastic process easily applicable in real time and on a long-term evaluation of azimuthal dependency of post-fit residuals under severe weather conditions. Data availability. GNSS data from the EUREF Permanent Network (EPN) stations are freely available through the anonymous FTP, e.g. from the EPN historical data centre at ftp://epncb.oma.be/ pub/obs/ maintained by the Royal Observatory of Belgium. Other GNSS and WVR data were primarily collected for the purpose of the COST Action ES1206 (GNSS4SWEC project; see Dousa et al. 2016) and cannot be distributed. Numerical weather model data fields from ALADIN-CZ were provided by the Czech Hydrometeorological Institute only for the purpose of the GNSS4SWEC project. Data were available only to the project members until the end of 2017 through the licence signed by the Research Institute of Geodesy, Topography and Cartography (RIGTC). The ERA-Interim data from the European Centre for Medium-Range Weather Forecasts (ECMWF, http://www.ecmwf.int/) were provided to GFZ by the German Weather Service. The Global Forecast System data were provided by the National Centers for Environmental Prediction (http://nomads.ncdc.noaa.gov/data/gfsanl). All the validation results in the form of figures and tables for all types of presented comparisons and stations can be provided by request to [email protected]. Competing interests. The authors declare that they have no conflict of interest. Acknowledgements. This study has been organized within the EU COST action ES1206 (GNSS4SWEC). The authors thank all the institutions that provided data for the benchmark campaign on which the validation was based on. Namely we want to thank S. Heise (GFZ) for providing the WVR data. The GFS data were provided by the National Centers for Environmental Prediction (www.ncep.noaa.gov). The ERA-Interim data were provided by the European Centre for Medium-Range Weather Forecasts (http://www.ecmwf.int/en/forecasts/datasets). M. Kaˇ cmaˇ rík, J. Douša and P. Václavovic acknowledge the support from the Czech Ministry of Education, Youth and Sports (project nos. LD14102 and LM2015079). E. Pottiaux (ROB) and H. Brenot (BIRA) acknowledge the support from the Solar-Terrestrial Centre of Excellence (STCE). J. Kapłon and P. Hordyniec (WUELS) acknowledge the support of Polish National Science Centre (project no. UMO-2013/11/D/ST10/03473) for financial support and Wrocław Center of Networking and Supercomputing (http://www.wcss.wroc.pl) for computational grant using MATLAB software license no. 101979. Edited by: J. Jones Reviewed by: two anonymous referees References Alber, C., Ware, R., Rocken, C., and Braun, J.: Obtaining single path phase delays from GPS double differences, Geophys. Res. Lett., 27, 2661–2664, https://doi.org/10.1029/2000gl011525, 2000. Askne, J. and Nordius, H.: Estimation of Tropospheric Delay for Microwaves from Surface Weather Data, Radio Sci., 22, 379– 386, https://doi.org/10.1029/rs022i003p00379, 1987. US Standard Atmosphere: US Standard Atmosphere, 1976, edited by: the National Oceanic and Atmospheric Administration, the National Aeronautics and Space Administration, and the US Air Force, US Government Printing Office, Washington, D. C., available at: http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/ 19770009539.pdf (last access: 8 June 2017), 1976. Bauer, H.-S., Wulfmeyer, V., Schwitalla, T., Zus, F., and Grzeschik, M.: Operational assimilation of GPS slant path delay measurements into the MM5 4DVAR system, Tellus A, 63, 263–282, https://doi.org/10.1111/j.1600-0870.2010.00489.x, 2011. Bender, M., Dick, G., Wickert, J., Schmidt, T., Shong, S., Gendt, G., Ge, M., and Rothacher, M.: Validation of GPS slant delays using water vapour radiometers and weather models, Meteorol. Z., 6, 807–812, https://doi.org/10.1127/0941-2948/2008/0341, 2008. Bender, M., Dick, G., Wickert, J., Ramatschi, M., Ge, M., Gendt, G., Rothacher, M., Raabe, A., and Tetzlaff, G.: Estimates of the information provided by GPS slant data observed in Germany regarding tomographic applications, J. Geophys. Res., 114, D06303, https://doi.org/10.1029/2008JD011008, 2009. Bender, M., Stosius, R., Zus, F., Dick, G., Wickert, J., and Raabe, A.: GNSS water vapour tomography: Expected improvements by combining GPS, GLONASS and Galileo observations, Adv. Space Res., 47, 886–897, https://doi.org/10.1016/j.asr.2010.09.011, 2011. Bender, M., Stephan, K., Schraff, C., Reich, H., Rhodin, A., and Potthast, R.: GPS Slant Delay Assimilation for Convective Scale NWP, Fifth International Symposium on Data Assimilation (ISDA), University of Reading, UK, July 18–22, 2016. Bennitt, E. and Jupp, A.: Operational Assimilation of GPS Zenith Total Delay Observations into the Met Office Numerical Weather Prediction Models, Mon. Weather Rev., 140, 2706– 2719, https://doi.org/10.1175/MWR-D-11-00156.1, 2012. Bevis, M., Businger, S., Chiswell, S., Herring, T., Anthes, R., Rocken, C., and Ware, R.: GPS Meteorology: Mapping Zenith Wet Delays onto Precipitable Water, J. Appl. Meteorol., 33, 379–386, https://doi.org/10.1175/15200450(1994)033<0379:gmmzwd>2.0.co;2, 1994. Böhm, J. and Schuh, H.: Vienna mapping functions, in: Proc. 16th Working Meeting on European VLBI for Geodesy and Astrometry, Leipzig, Germany, Verlag des Bundesamtes für Kartographie und Geodäsie, 131–143, 2003. Böhm, J., Werl, B., and Schuh, H.: Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range Weather Forecasts opAtmos. Meas. Tech., 10, 2183–2208, 2017 www.atmos-meas-tech.net/10/2183/2017/ M. Kaˇ cmaˇ rík et al.: Inter-technique validation of tropospheric slant total delays 2207 erational analysis data, J. Geophys. Res., 111, B02406, https://doi.org/10.1029/2005JB003629, 2006a. Böhm, J., Niell, A., Tregoning, P., and Schuh, H.: Global MappingFunction (GMF): A new empirical mapping function based on numerical weather model data, Geophys. Res. Lett., 33, L07304, https://doi.org/10.1029/2005GL025546, 2006b. Braun, J., Rocken, C., and Ware, R.: Validation of line-of-sight water vapour measurements with GPS, Radio Sci., 36, 459–472, https://doi.org/10.1029/2000rs002353, 2001. Braun, J., Rocken, C., and Liljergen, J.: Comparisons of Lineof-Sight Water Vapor Observations Using the Global Positioning System and a Pointing Microwave Radiometer, J. Atmos. Ocean. Tech., 20, 606–612, https://doi.org/10.1175/15200426(2003)20<606:colosw>2.0.co;2, 2002. Brenot, H.: Potential of ground-based GPS measurements for the study of Mediterranean heavy rains, PhD of the French State, http://tel.archives-ouvertes.fr/tel-00012085 (last access: 7 June 2017), 2006. Brenot, H., Neméghaire, J., Delobbe, L., Clerbaux, N., De Meutter, P., Deckmyn, A., Delcloo, A., Frappez, L., and Van Roozendael, M.: Preliminary signs of the initiation of deep convection by GNSS, Atmos. Chem. Phys., 13, 5425–5449, https://doi.org/10.5194/acp-13-5425-2013, 2013. Chen, G. and Herring, T. A.: Effects of atmospheric azimuthal asymmetry on the analysis of space geodetic data, J. Geophys. Res., 102, 20489–20502, https://doi.org/10.1029/97JB01739, 1997. Dach, R., Lutz, S., Walser, P., and Fridez, P. (Eds.): Bernese GNSS Software Version 5.2. User manual, Astronomical Institute, University of Bern, Bern Open Publishing, https://doi.org/10.7892/boris.72297, 2015. Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: Configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597, https://doi.org/10.1002/qj.828, 2011. Deng, Z., Bender, M., Zus, F., Ge, M., Dick, G., Ramatschi, M., Wickert, J., Löhnert, U., and Schön, S.: Validation of tropospheric slant path delays derived from single and dual frequency GPS receivers, Radio Sci., 46, RS6007, https://doi.org/10.1029/2011RS004687, 2011. Douša, J. and Václavovic, P.: Real-time zenith tropospheric delays in support of numerical weather prediction applications. Adv. Space Res., 53, 1347–1358, https://doi.org/10.1016/j.asr.2014.02.021, 2014. Douša, J., Dick, G., Kacmarík, M., Brožková, R., Zus, F., Brenot, H., Stoycheva, A., Möller, G., and Kaplon, J.: Benchmark campaign and case study episode in central Europe for development and assessment of advanced GNSS tropospheric models and products, Atmos. Meas. Tech., 9, 2989–3008, https://doi.org/10.5194/amt-9-2989-2016, 2016. De Haan, S., Marel, van der H., and Barlag, S.: Comparison of GPS slant delay measurements to a numerical model: case study of a cold front passage, Phys. Chem. Earth, 27, 317–322, https://doi.org/10.1016/s1474-7065(02)00006-2, 2002. Flores, A., Rius, A., Vilá-Guearou, J., and Escudero, A.: Spatio-temporal tomography of the lower troposphere using GPS signals, Phys. Chem. Earth A, 26, 405–411, https://doi.org/10.1016/s1464-1895(01)00074-6, 2001. Gradinarsky, L. P.: Sensing Atmospheric Water Vapor Using Radio Waves, PhD thesis, School of Electrical Engineering, Chalmers University of Technology, Göteborg, Sweden, 2002. Guerova, G., Bettems, J. M., Brockmann, E., and Matzler, C.: Assimilation of COST 716 Near-Real Time GPS data in the nonhydrostatic limited area model used at MeteoSwiss, Meteorol. Atmos. Phys., 91, 149–164, https://doi.org/10.1007/s00703-0050110-6, 2006. Guerova, G., Jones, J., Douša, J., Dick, G., de Haan, S., Pottiaux, E., Bock, O., Pacione, R., Elgered, G., Vedel, H., and Bender, M.: Review of the state of the art and future prospects of the groundbased GNSS meteorology in Europe, Atmos. Meas. Tech., 9, 5385–5406, https://doi.org/10.5194/amt-9-5385-2016, 2016. Ha, S.-Y, Kuo, Y.-H., Guo, Y.-R., Rocken, C., and Van Hove, T.: Comparison of GPS slant wet delay measurements with model simulations during the passage of a squall line, Geophys. Res. Lett., 29, https://doi.org/10.1029/2002gl015891, 2002. Hobiger, T., Ichikawa, R., Koyama, Y., and Kondo, T.: Fast and accurate ray-tracing algorithms for real-time space geodetic applications using numerical weather models, J. Geophys. Res.- Atmos., 113, D20302, https://doi.org/10.1029/2008jd010503, 2008. Kaˇ cmaˇ rík, M., Douša, J., and Zapletal, J.: Comparison of GPS slant wet delays acquired by different techniques, Acta Geodyn. Geomater., 9, 427–433, 2012. Kawabata, T., Shoji, Y., Seko, H., and Saito, K.: A Numerical Study on a Mesoscale Convective System over a Subtropical Island with 4D-Var Assimilation of GPS Slant Total Delays, J. Meteorol. Soc. Jpn., 91, 705–721, https://doi.org/10.2151/jmsj.2013510, 2013. Lemoine, F. G., Kenyon, S. C., Factor, J. K., Trimmer, R. G., Pavlis, N. K., Chinn, D. S., Cox, C. M., Klosko, S. M., Luthcke, S. B., Torrence, M. H., Wang, Y. M., Williamson, R. G., Pavlis, E. C., Rapp, R. H., and Olson T. R.: The development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96, NASA/TP-1998-206861, https://doi.org/10.1007/978-3-662-03482-8_62, 1998. Li, X., Zus, F., Lu, C., Dick, G., Ning, T., Ge, M., Wickert, J., and Schuh, H.: Retrieving of atmospheric parameters from multi-GNSS in real time: Validation with water vapor radiometer and numerical weather model, J. Geophys. Res.-Atmos., 120, https://doi.org/10.1002/2015JD023454, 2015a. Li, X., Zus, F., Lu, C., Ning, T., Dick, G., Ge, M., Wickert, J., and Schuh, H.: Retrieving high-resolution tropospheric gradients from multiconstellation GNSS observations. Geophys. Res. Lett., 42, 4173–4181, https://doi.org/10.1002/2015GL063856, 2015b. MacMillan, D. S.: Atmospheric gradients from very long baseline interferometry observations, Geophys. Res. Lett., 22, 1041– 1044, https://doi.org/10.1029/95GL00887, 1995. Mahfouf, J.-F., Ahmed, F., Moll, P., and Teferle, F. N.: Assimilation of zenith total delays in the AROME France convecwww.atmos-meas-tech.net/10/2183/2017/ Atmos. Meas. Tech., 10, 2183–2208, 2017