Atmospheric Research 138 (2014) 414–426 Contents lists available at ScienceDirect Atmospheric Research journal homepage: www.elsevier.com/locate/atmos Ground-based GNSS ZTD/IWV estimation system for numerical weather prediction in challenging weather conditions Witold Rohm a,b,⁎, Yubin Yuan a, Bertukan Biadeglgne c, Kefei Zhang a, John Le Marshall c a b c RMIT University, SPACE Research Centre, GPO Box 2476 Melbourne, VIC 3001, Australia Wroclaw University of Environmental and Life Sciences, Institute of Geodesy and Geoinformatics, Grunwaldzka 53, 50-357 Wroclaw, Poland Australian Bureau of Meteorology, Melbourne, VIC, Australia a r t i c l e i n f o Article history: Received 2 August 2013 Received in revised form 23 October 2013 Accepted 30 November 2013 Keywords: ZTD/IWV estimation GNSS meteorology Verification a b s t r a c t The Global Navigation Satellite Systems (GNSS) are one of the very few tools that can provide continuous, unbiased, precise and robust atmosphere condition information. The extensive research of GNSS space-based segment (e.g. available precise, real-time satellite orbits and clocks), unlimited access to the ground-based Continuously Operating Reference Stations (CORS) GNSS networks along with the well established data processing methods provides an unprecedented opportunity to study the environmental impacts on the GNSS signal propagation. GNSS measurements have been successfully used in precise positioning, tectonic plate monitoring, ionosphere studies and troposphere monitoring. However all GNSS signals recorded on the ground by CORS are subject to ionosphere delay, troposphere delay, multipath and signal strength loss. Nowadays, the GNSS signal delays are gradually incorporated into the numerical weather prediction (NWP) models. Usually the Zenith Total Delay (ZTD) or Integrated Water Vapour (IWV) have been considered as an important source of water vapour contents and assimilated into the NWP models. However, successful assimilation of these products requires strict accuracy assessment, especially in the challenging severe weather conditions. In this study a number of GNSS signal processing strategies have been verified to obtain the best possible estimates of troposphere delays using a selection of International GNSS Service (IGS) orbit and clock products. Three different severe weather events (severe storm, flash flooding, flooding) have been investigated in this paper. The strategies considered are; 1) Double Differenced (DD) network solution with shortest baselines, 2) DD network solution with longest baselines, 3) DD baseline-by-baseline solution (tested but not considered), 4) Zero Differenced (ZD) Precise Point Positioning (PPP) based on ambiguity float solutions, all with precise orbits and clocks, and real time clocks and predicted orbits. The quality of the estimates obtained has been evaluated against radiosonde measurements, Automatic Weather Station (AWS) observations, NWP (assimilation step without ground-based GNSS data) and ZTD estimates from the well established IGS processing centre, the Center of Orbit Determination in Europe (CODE). It shows that the ZTD and IWV estimates from the DD short baseline solution are robust with usually a very small bias (−2.7 to −0.8 mm) and errors of less than 10 mm (7.6–8.5) (ZTD) or 3 mm (2.6–2.7) (IWV). The DD short baseline network solution was found to be the most reliable method in the considered case studies, regardless of the type of orbits and clocks applied. © 2013 Elsevier B.V. All rights reserved. 1. Introduction ⁎ Corresponding author at: Wroclaw University of Environmental and Life Sciences, Institute of Geodesy and Geoinformatics, Grunwaldzka 53, 50-357 Wroclaw, Poland. E-mail address: witold.rohm@up.wroc.pl (W. Rohm). 0169-8095/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.atmosres.2013.11.026 The successful assimilation of GNSS space-based observations from Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) mission and other Low W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 Earth Orbiters (LEO) satellites in NWP models (Marshall et al., 2010) has significantly stimulated the interest of meteorological community in using ground-based GNSS observations to improve weather forecasting. The ground-based GNSS observations are collected by dense networks of CORS and processed by a central unit to obtain precise atmosphere information to improve the quality of satellite clocks and orbits and therefore improve positioning solutions. The by-product of the CORS data processing is the troposphere delay linked with the weather conditions (pressure, temperature and water vapour). The most important parameter in severe weather forecasting is the water vapour (WV) content since it is a WV induced phenomenon. The WV plays an active role in energy exchange between climatic zones, and in a finer scale as a vertical channel of energy transfer (e.g. by latent heat). Therefore, all major severe weather events that are linked to the WV can be detected by GNSS signals. However, the rapidly changing weather conditions such as those presented before, during and after severe weather events pose a big challenge for a successful data processing of GNSS signals to achieve the required accuracy of troposphere conditions. The ZTD accuracy requirements imposed by NWP assimilation strategy as reported by (Barlag et al., 2004) are 10 mm, whereas the IWV should be delivered with an accuracy better than 5 mm. This study is focused on identifying potential problems of GNSS data processing in challenging weather conditions, quantifying the magnitude of errors introduced by such phenomena and eventually removing their impacts. The ultimate goal of this study is to provide optimal strategy to retrieve troposphere phase delay from GNSS ground based observations for NWP assimilation algorithms. The GNSS signal is bended, attenuated and delayed in the atmosphere, on two basic layers through which it propagates predominately, troposphere and ionosphere. The signal bending effect, widely used in satellite-to-satellite GNSS data processing, has negligible impact on the GNSS signals recorded by the ground-based receivers. The neutral refractivity N0, a main cause of delay, could be split into the following form (Thayer, 1974; Solheim et al., 1999): N0 ¼ N d þ Nv þ Nnon þ Ndiff ð1Þ where Nd is the dry delay and Nv is the water vapour delay, whereas Nnon and Ndiff are linked with the non-gaseous parts of atmosphere (like dust) and hydrometeors, respectively. The last two elements of Eq. (1) have usually very limited impacts on refractivity and in practice they are not considered in data processing strategies. The neutral refractivity is calculated from the equation below: N0 ¼ Z d −1
k1
 pd e e −1 þ Zv k2 þ k3 2 ; T T T 415 The dynamics of the changes in WV content e and temperature T will result in substantial time and space variations of the wet refractivity Nv in time and space (Eq. (2)). In severe weather conditions the changes in the wet refractivity are significant (Manning et al., 2012). Consequently, there are several important issues linked with GNSS data processing for troposphere: 1) the azimuthal inhomogeneity of the troposphere, 2) the unknown and variable correlation time of the troposphere conditions, and 3) a high value of stochastic parameter reflecting the variation of the troposphere conditions. The propagation environment for GNSS signals arriving from different directions (both elevation angles and azimuths) varies, hence the functional model of phase propagation delay should reflect these inhomogeneities. Usually to resolve this issue troposphere gradients are introduced. Another important factor related to the dynamics of the troposphere is the troposphere delay parametrisation time step size. The linear parametrisation (Dach et al., 2007) assumes that the path delay is constant within the time step within the specified a priori standard deviation. Therefore the stochastic modelling of signal phase delay should be carefully considered, because the a posteriori ZTD error will reflect how the functional model fits into the data. The GNSS signal propagation from a satellite to the receiver through the neutral atmosphere is subject to the change of propagation speed expressed as a Slant Total Delay (STD) (δρR(zRS)) and given by the following equation:     −6 −6 S δρR zR ¼ 10
∫N0 ds ¼ 10
∫Nd ds þ ∫Nv ds : ð3Þ A good approximation of STD can be expressed usually, as a first a priori value (δρapr,R), according to the following formula:       S S S δρapr;R zR ¼ mapH zR
ZHD þ mapW zR
ZWD; ð4Þ where mapH(zRS) is the mapping function for the dry delay ZHD and mapW(zRS) is the mapping function for the wet delay ZWD. The STD is not related to the frequency of the GNSS signals, therefore the impact of troposphere, unlike the ionosphere, cannot be reduced by linear combinations of the signal frequencies. Usually to derive signal's delay in the neutral atmosphere (Dach et al., 2007; Herring et al., 2010), corrections at stations to the a priori model along with station coordinates are estimated. The functional model adopted is as follows (Dach et al., 2007):     ∂mapN S S h S n S δρR ¼ δρapr;R zR þ δ ρR ðt Þ
mapN zR þ δ ρR ðt Þ cosAR þ …… ∂z ∂mapN S e sinAR ; þ δ ρR ðt Þ ∂z ð5Þ ð2Þ where pd is a dry pressure (the pressure of atmosphere excluding water vapour partial pressure) and temperature T, water vapour partial pressure e coupled with temperature T. −1 −1 Z d is an inversion of dry compressibility factor and Z v is an inversion of wet compressibility factor (Owens, 1967). The k1, k2, and k3 parameters are atmospheric refractivity constants given by (Kleijer, 2004). where, zRS is the zenith angle of satellite S as seen from station R, ASR is a satellite azimuth, δhρR(t) is a time t dependent delay in zenith direction at a point R, mapN(zRS) is a mapping function, δnρR(t) is a time dependent gradient in north direction, and δeρR(t) is a time dependent gradient in east direction. However, the realization of this principle varies between data processing software packages and processing strategies. The stochastic of the troposphere parameter is rather simple (Schüler, 2001). The absolute value of δρR can be 416 W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 constrained to an arbitrary value (absolute constraining), or alternatively the rate of change of troposphere delay at the station δρδt could be kept fixed. The second option is to use a maximum correlation time. These stochastic parameters are linked with troposphere conditions, therefore their empirical evaluation is presented in the case study section. The ZWD, as a measure of WV content and temperature conditions in the troposphere, is derived from the following equation: R h ZWD ¼ δ ρR ðt Þ−ZHD: ð6Þ The relation between ZWD and the WV content in the atmosphere is expressed by IWV and given by the following equation (i.e. (Kleijer, 2004)): IWV ¼ ZWD

′ k2 þ k3 TM −1 10−6
Rw ¼ ZWD
Q ð7Þ where Rw = 461.525 ± 0.003 [J kg−1 K−1] is the specific gas constant for water vapour, k2′ = 24 ± 11 [K hPa−1], and k3 = 3.75 ± 0.03 [105 K2 hPa−1] are refraction constants (Schüler, 2001) and TM ≈ 70.2 ± 0.72 ⋅ T0 is the weighted mean water vapour temperature of the atmosphere, and T0 is the surface temperature (Mendes and Langley, 1999). Alternatively TM value can be retrieved from the profile of temperature T and water vapour partial pressure e according to the equation: TM ¼ ∑Te : ∑Te ð8Þ 2 This study focuses on the investigation of the most robust yet efficient method to derive ZTD and IWV in severe weather conditions. The challenging weather conditions will presumably result in the following propagation effects: 1. Non-homogeneity of the horizontal distribution of the atmosphere masses. To account for this effect, the ZTD gradients are introduced in Eq. (5). The variation of the wet refractivity in the approach demonstrated in this study is measured with a time autocorrelation function of wet refractivity distribution in selected troposphere altitude layers, in normal and storm prone conditions. 2. The increased ZTD estimation uncertainty as a response to rapid changes in temperature and WV content and the parametrisation scheme with constant ZTD for the epoch to be estimated. This effect is investigated by comparing the estimation uncertainty of ZTD in normal conditions to ZTD in the severe storm conditions. 3. Low observation signal-to-noise ratio during the storm activity. The impact of atmospheric loss fluctuations on the estimated parameters is challenging to quantify (Misra and Enge, 2001; Christian et al., 2004) and has not been discussed in this study. 4. In DD strategies, uneven weather conditions on both sides of the baseline increase the reliability of ZTD absolute estimation. This effect is discussed in the following sections. In order to select the most appropriate approach, several processing strategies (in Section 2) with final and real-time IGS products (clocks and orbits) have been validated. In DD processing mode, three data processing scenarios are examined, these include: 1) the network solution with shortest baselines, 2) the network solution with longest baselines, and 3) baseline-by-baseline estimation with shortest baselines. The ZD phase observations (PPP) have also been investigated in this study. Majority of the computing work has been performed using Bernese GPS Software version 5.0. However the methodology discussed in this manuscript is equally applicable to any GNSS data processing software package. The data used for control and validation (contained in Section 3) of results comprise: deterministic models Global Pressure and Temperature (GPT) (Boehm et al., 2007), UNB3m (Leandro et al., 2008), ground-based Automatic Weather Stations managed by Australian Bureau of Meteorology (BoM), radiosonde observations from 14 locations in Australia and Australian Community Climate and Earth-System Simulator —Regional (ACCESS-R) outputs of analysis step. Section 4 contains the analysis of three severe case studies. The case studies are composed of: 2010 severe storm that hit the Melbourne metropolitan area, 2011 flooding in Victoria and 2012 flash flooding in Victoria. The paper is concluded in Section 5. 2. ZTD estimation The standard procedure to obtain high quality coordinates requires double differencing of phase observations to remove the satellite and receiver clock errors and exploit the integer nature of ambiguities (Hofmann-Wellenhof et al., 2008). Usually the final solution is obtained on Ionosphere Free combination (L3) (Dach et al., 2007). Procedure adopted in this study in case of DD processing follows (Bosy et al., 2003) with minor modifications. Precedence to the ZTD estimation, due to a large number of stations (e.g. 218 in 2012), the normal equation set-up is separated into 3 regions. Within each region the same procedure is applied as follows: • • • • • • • • the reference station coordinates are propagated into observation epoch, the Lagrange polynomials are fit into the discrete satellite positions, the receiver clocks are synchronised using code-based standard positioning solution, the single differencing of the phase observations is performed using predefined baselines (skeleton network only) and selected differentiation strategy (shortest, longest), the phase observations undergo quality control including: cycle slip detection, marking and removing spurious observations, the first estimation of float ambiguities on wide —lane (L5) wavelength is performed, the second estimation of integer values of ambiguities on narrow —lane (L3) wavelength is obtained, the previously estimated ambiguities are introduced to form and store normal equations. W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 417 This procedure is followed by steps common for the whole network: • combination of daily normal equation files across the whole network into one file per network per day • estimation of weekly coordinates using constrained translation parameters based on velocities and coordinates of IGS reference stations, • estimation of troposphere parameters (ZTD, and troposphere gradients) with daily normal equation files with minimum constraint conditions imposed on IGS stations' coordinate translation parameters. The last item listed above contains several standardised steps and the procedure given below is equally valid for Zero Differenced (ZD) observations. According to Eq. (2) refractivity is in 80–90% linked with pressure of dry atmosphere pd, and in 10–20 % to the distribution of water vapour e and temperature T and consequently to the stations height, so the first step would be the determination of the total delay applying globally adopted mean of meteorological parameters with corrections to the stations height. Then, at each station ZTD with gradient parameters is modelled as a correction to the a priori model given previously, according to Eq. (5). The IGS Final orbits and clocks (PREC) are combined from a number of contributing IGS Analysis Centres, using six, largely independent, software packages (i.e. BERNESE, GAMIT, GIPSY, NAPEOS, EPOS and PAGES). The IGS Final orbit/clocks are usually available in two weeks after the last observation. The anticipated accuracies of final orbits and clock products are 1–2 cm and 0.02–0.06 ns respectively (Kouba, 2009). The real-time clocks and orbits (PRED) are the outcomes of the IGS Real-time Pilot Project (IGS-RTPP). The IGS real-time processing centre estimates and distributes high-resolution GPS corrections of broadcast orbits and high frequency clocks for PPP applications. The orbits used in this study are ultrarapid predictions from 8 processing centres, the anticipated accuracy of these orbits is 5 cm. The clock product is estimated in the real time mode (Li et al., 2013) from the network of global IGS stations, the quality achieved is on the order of 0.10 ns (http://rts.igs.org/monitor/). 2.1. Double Difference processing The standard equation for DD phase observations reads as follows (Hofmann-Wellenhof et al., 2008): ST ST ST ST ST λ
ΦPR ðt Þ ¼ δr PR þ λ
NPR þ δρPR þ λ
PR : ð9Þ The final ZTD estimation step for DD case is conducted according to Dach et al. (2007), Kubo et al. (2012), and Strang and Borre (1997). The linearised Eq. (9) in a matrix form read: AX = l where A represents the design matrix, X is the parameter to be estimated and l is a set of observations. The estimation of troposphere parameters is obtained with previously derived, integer ambiguities N, therefore in Eqs. (10)–(12) N values do not appear. The following is an example of two receivers, four satellites for a single epoch in DD case. 2 S1 T 2 6 aX PR ðt Þ 6 6 A ¼ 6 aSX1 T 3 ðt Þ 6 PR 4 S T aX1PR 4 ðt Þ S T aZ1PR 2 ðt Þ S T aZ1PR 3 ðt Þ S T aZ1PR 4 ðt Þ aY1PR 2 ðt Þ aY1PR 3 ðt Þ aY1PR 4 ðt Þ S T S T S T   ∂mapN S T S T δmapN zPR1 2 δ cosAPR1 2 ∂z   ∂mapN S T S T δmapN zPR1 3 δ cosAPR1 3 ∂z   ∂mapN S T S T δ δmapN zPR1 4 cosAPR1 4 ∂z 3 ∂mapN S T sinAPR1 2 7 ∂z 7 ∂mapN S T 7 δ sinAPR1 3 7 7 ∂z 5 ∂mapN S1 T 4 δ sinAPR ∂z δ ð10Þ 3 Δ X PR 6 ΔY PR 7 7 6 7 6 ΔZ PR 7 6 7 6 h 6 δ ρP ðt Þ 7 7 6 n X ¼ 6 δ ρP ðt Þ 7 7 6 e 6 δ ρP ðt Þ 7 7 6 h 6 δ ρR ðt Þ 7 7 6 n 4 δ ρR ðt Þ 5 e δ ρR ðt Þ ð11Þ  3 S T S T S T λ
ΦPR1 2 ðt Þ−δr PR1 2 −Δδρapr;R zPR1 2 6  7 S T S T S T 7 6 l ¼ 6 λ
ΦPR1 3 ðt Þ−δr PR1 3 −Δδρapr;R zPR1 3 7: 4  5 S1 T 4 S1 T 4 S1 T 4 λ
ΦPR ðt Þ−δr PR −Δδρapr;R zPR ð12Þ 2 2 The troposphere parameter is estimated as a piecewise linear function with monotonic parameters in zenith direction δhρR(t) along with gradient parameters δnρR(t) and δeρR(t). Several, more specific strategies were tested in the processing case studies and are listed below. 2.1.1. Baseline-by-baseline solution The investigation of baseline length impact has been performed using baseline-by-baseline processing. In that case, for each shortest baseline, as given by Bernese GPS Software (Dach et al., 2007), design matrix A has been setup for whole day, covering 418 W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 satellites and two stations. The shortest baseline strategy produces vectors of the mean distance of 70 km. The ZTDs' distance dependent correlation coefficients (from PPP solution), for baselines shorter than 100 km in flat lowlands (station MOBS 41 m), are above 0.95. If the highest point in the network would be considered (station HOTH 1774 m), then for the same interstation distance correlation coefficient would be above 0.90. Hence there are similar troposphere conditions on both ends of the baseline. Therefore the last 3 columns of matrix A (10) are close to 0, which produces a high uncertainty of solutions and amplification of errors presented in the observations. The mean discrepancy between reference data and outcomes of this strategy is 9 mm and 31 mm in terms of bias and standard deviation, respectively. The maximum outliers reach up to 400 mm. Therefore baseline-by-baseline solution has not been regarded as an appropriate strategy for troposphere estimation and is not presented in our figures or tables. 2.1.2. Shortest baselines network solution The other approach, routinely adopted to derive precise coordinates and velocities of national and regional reference networks (Bosy et al., 2009), requires setting up a network of n − 1 independent shortest baselines (n is the number of points) (Dach et al., 2007). The matrix A (10) previously used to obtain troposphere solutions for one baseline (PR) is converted to a stacked matrix As covering the whole network. Observation matrix ls is composed of all observations in the time span considered (usually a day) across all baselines. The vector of unknown Xs consists of coordinates, zenith delays and gradient parameters in the north and east directions. Therefore, the number of unknowns is limited to: 3 ⋅ nst number of baselines in the network (one for each coordinate), plus nztd times the interval length of ZTD estimation plus one (extra term in the beginning or in the end of the observation window) and ngrad times interval length of the ZTD gradient estimation (usually much lower then nztd). The key issue of all Weighted Least Squares (WLS) (Strang and Borre, 1997) estimation is covariance matrix C. The mathematical correlation between DD observations (Schüler, 2006) has been considered (Dach et al., 2007), as well as correlations between baselines (Hofmann-Wellenhof et al., 2008). In addition, minimum constraints have been applied on translation parameters and an elevation  dependent weighting σ0() function has been applied. The campaigns discussed in this study have been processed by applying ZTD relative constraints (σrelZTD) with maximum time constraints. The values associated with relative constraints are set based on the empirical tests using ZTD autocovariance. In the Bernese GPS software, constraints are introduced as fictitious observations (Dach et al., 2007) yfictitious = H ⋅ x where H is a matrix of ones or zeros when an observation add 1 or not 0. All statistical considerations given above translate into variance covariance matrix Cs of the following form (Hofmann-Wellenhof et al., 2008; Schüler, 2006): 4 62 6 62 ¼ σ 0 ðÞ
6 62 6 41 1 2 C DD 2 4 2 1 2 1 2 2 4 1 1 2 2 1 1 4 2 2 1 2 1 2 4 2 3 1 17 7 27 7 27 7 25 4 ð13Þ Fig. 1. Longest baseline selection scheme. Boxes in red represent the selected stations in set A; boxes in green represent the selected stations in set B; boxes in blue represent the stations to be selected in set B. W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 C DD 6 0 6 6 0 6 6 0 CS ¼ 6 6 0 6 6 0 6 4 0 0 2 0 σX 0 0 0 0 0 0 0 0 σY 0 0 0 0 0 0 0 0 σZ 0 0 0 0 0 0 0 0 σ relZTD 0 0 0 0 0 0 0 0 σ relZTD 0 0 0 0 0 0 0 0 σ δn 0 3 0 0 7 7 0 7 7 0 7 7: 0 7 7 0 7 7 0 5 σ δe 419 ð14Þ The covariance matrix CS is inverted to obtain the weighting matrix P. The DD network solution discussed in this paragraph considers full correlation matrix (Eq. (14)) between baselines and within double difference observations (Eq. (13)). This extensive correlation modelling and the CORS network size, spanning thousands of kilometres from Adelaide (South Australia) to Townswille (North Queensland) lead to absolute delay estimation, regardless of the baseline length. 2.1.3. Longest baselines network solution The longest baseline solutions are also evaluated to investigate the best solution for GNSS data processing. In this study the testing approach used for the longest baselines is demonstrated in Fig. 1 using the following steps: Step 1 All the n stations are located in set B while the selected stations consist of set A which is originally empty (with zero value elements); Step 2 We start from selecting one random station A[0]. Distances between A[0] and the other n − 1 stations are compared. The station with the longest baseline is where A[0] will be selected next. After this step there are A[0] and A[1] in set A while the rest n − 2 stations are still in set B; Step 3 Similar to Step 2, the length of baselines between A[1] and the stations in set B is compared. This step is repeated iteratively until no stations are left in set B; Step 4 The set B is re-ordered as set A in which two adjacent stations form a baseline. This strategy creates the baselines of 900 km average length. In this strategy the baseline lengths are significantly longer (usually 10 times) than those in the network shortest solution. However the mathematical observation models remain unchanged (10, 11, 12). The most significant differences are threefold: 1) the last columns of matrix A (10) contains coefficients that are uncorrelated (distance approx. 500 km); 2) as a consequence, the estimation of ambiguities is more challenging, due to the shortest common observation window; and 3) different ionosphere conditions on both ends of the baseline. Therefore, the number of the longest baseline solution benefits (no correlation of troposphere conditions) are outnumbered by disadvantages linked with station separation (ambiguities fixing problem, lower number of observations). The network solution regardless of the configuration (shortest or longest) shows actually the same performance in terms of estimating troposphere conditions, which has been proved in the case study section. 2.2. Precise Point Positioning processing The last method considered in this study is based on PPP processing, of ZD phase observations. The flow is similar to DD processing, already discussed in great detail in the previous section. In the Bernese GPS Software the ambiguities of ZD observations are not fixed to the integer values. Studies show (Sunil Bisnath —personal communication, this paper) that PPP float ambiguities processing of long GNSS observation session will have the same accuracy of troposphere estimates as the one with fixed ambiguities. Similar to the DD processing, the float ambiguity N are estimated and removed from processing priori to the ZTD estimation, therefore N values are omitted in Eqs. (16)–(18). The ZD phase observation model is slightly different to the DD one (Eq. (9)), the commonly adopted model follows (Dach et al., 2007) and it reads as follows: S S S S λ
ΦR ðt Þ ¼ δr R þ λ
N R þ δρR þ c
R ðt Þ: ð15Þ The comparison of the PPP fundamental Eq. (15) and the DD basic Eq. (9) reveals an extra term c ⋅ R(t) that is linked with the receiver clock error (Eq. (15)). The satellite clock bias (PREC and PRED product) is used to correct the phase observations (Eq. (15)). The time series of this parameter exhibits short-term variations super imposed over a long-term trend (Collins et al., 2010). It is also important to note that all quantities given in Eq. (15), unlike in DD model (Eq. (9)) (including ambiguities) are absolute values. Therefore any biases presented in the orbits (Douša, 2010) and satellite clocks will have strong signature in the 420 W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 Fig. 2. The scheme presenting the ZTD/IWV evaluation process using AWS, RS and NWP. troposphere values obtained. Eq. (15) is linearised and (in the case of four satellites in one epoch) it reads as follows l = A ⋅ X where: 2 1 6 aX R ðt Þ 6 6 2 6 aX ðt Þ 6 R A¼6 6 3 6 aX R ðt Þ 6 4 4 aX R ðt Þ 1 aY R ðt Þ 2 aY R ðt Þ 3 aY R ðt Þ 4 aY R ðt Þ   1 S aZ R ðt Þ 1 mapN zR1   2 S aZ R ðt Þ 1 mapN zR2   3 S aZ R ðt Þ 1 mapN zR3   4 S aZ R ðt Þ 1 mapN zR4 ∂mapN S cos AR1 ∂z ∂mapN S cos AR2 ∂z ∂mapN S cos AR3 ∂z ∂mapN S cos AR4 ∂z 3 ∂mapN S sinAR1 7 ∂z 7 ∂mapN S 7 sinAR2 7 7 ∂z 7 ∂mapN S 7 sinAR3 7 7 ∂z 5 ∂mapN S4 sinAR ∂z ð16Þ 3 ΔX R 6 ΔY R 7 7 6 6 ΔZ R 7 7 6 7 X¼6 6 c h
R ðt Þ 7 6 δ ρ ðt Þ 7 6 n S 7 4 δ ρ ðt Þ 5 ð17Þ  3 S S S λ
ΦR1 ðt Þ−δr R1 −δρapr;R zR1  7 6 S2 7 6 λ
ΦS2 ðt Þ−δr S2 −δρ 7 6 apr;R zR R R   7: l¼6 7 6 S3 S3 S 6 λ
ΦR ðt Þ−δr R −δρapr;R zR3 7 4  5 S S S λ
ΦR4 ðt Þ−δr R4 −δρapr;R zR4 ð18Þ 2 e S δ ρS ðt Þ 2 The essential difference between parameter vector X in PPP (Eq. (17)) and in DD (Eq. (17)), is that the number of estimated parameters in the case of PPP is larger by at least one parameter per epoch, which is the receiver clock term cR(t). Hence the design matrix A in PPP (Eq. (16)) has one column more than matrix A in DD (Eq. (10)) case. The unwanted correlation between columns of A will not occur, unlike in the DD case, because matrices A (Eq. (16)) and l (Eq. (18)) contains observations from independent satellites. 3. Validation of ZTD and IWV 3.1. GNSS reference data The data processed in this study were validated against reference GNSS troposphere and meteorological products, for both ZTD and IWV. The verification process follows the scheme presented in Fig. 2. The data processed by the Center for Orbit Determination in Europe (CODE) Astronomical Institute University of Bern (AIUB) has been used as a reference, mainly final troposphere 421 W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 solutions for Australian stations. The processing strategy applied is summarised as follows1: Step 1 The basic observable is GPS/GLONASS carrier phase; code is only used for receiver clock synchronisation, elevation angle cutoff is 3 degrees, sampling rate is 3 minutes, Step 2 continuous parameterisation of troposphere ZTD and horizontal gradient parameters, allowing for connection of the parameters at day boundaries, Step 3 a priori model for troposphere Saastamoinen-based hydrostatic (using GPT) mapped with the dry-GMF, Step 4 zenith delay corrections are estimated based on the wet Global Mapping Function (GMF) (Boehm et al., 2006) in 2 h intervals. N– and E–W horizontal delay parameters are solved for every 24 h, Step 5 both zenith and gradient parameters are treated as completely unconstrained, Step 6 ITRF2005/ITRF2008 reference frame for fiducial IGS stations and orbits, Step 7 Datum definition: 3 no-net translation conditions, 3 no-net rotation conditions, geocenter coordinates constrained nominally to zero values. 3.2. Meteorological data The set of meteorological sensors and models includes: • The AWS, operated by the Bureau of Meteorology, are designed to provide real-time data for weather forecasting and warning as well as high quality data for climate applications. AWS measure and report meteorological variables such as temperature TAWS, humidity RHAWS, pressure pAWS, wind and rainfall. Only the first three parameters are used in this study. The selection of useful stations was based on the distance to GNSS station. In this study a 50-km radius has been adopted. • The radiosonde (RS) vertical profile measurements from 17 stations across eastern Australia are provided every 12 or 24 h. However, only 14 stations have been used, due to the 100 km separation limit between RS and GNSS station. The vertical profiles comprise: pressure on the mandatory and significant levels, temperature, dew point temperature (humidity equivalent), height, wind speed and wind direction. Only the pressure, temperature and dew point temperature from mandatory levels have been utilised in this study. • The NWP analysis outputs from the Australian Community Climate and Earth System Simulator (ACCESS-A) model were used to validate the GNSS data. The NWP model outputs (analysis run) with the time resolution of 6 h (for base time: 00, 06, 12, 18 UTC) using a centred 6 hour observational data window. The model covers a spatial outline from − 45.0 ° S to 4.73 ° N and 95oE to 169.9 ° E. The ACCESS system is based on the UK Meteorological office Unified Model, and uses a number of data sources to produce forecasts (e.g., AIRS radiances, COSMIC, radiosondes, pilot and profiler winds, ATVOS radiances, scatterometer winds, AMDAR, land surface stations) (Marshall et al., 2010). The model in the horizontal plain 1 http://igscb.jpl.nasa.gov/igscb/center/analysis/archive/code_20080528. acn. contains 229 nodes with the grid spacing of 0.11 ° (12 km) and in the vertical direction it utilises terrain following hybrid (pressure/height) coordinates with 48 levels. The model's forecast and analysis covers a number of meteorological parameters such as zonal wind speed, meridonal wind speed, cloud cover, temperature, precipitation etc. This study however, considers only pressure, temperature and water vapour partial pressure (given as a mixing ratio). • The 3D wet refractivities, outputs of tomography modelling, for 6 days in March 2010, were used to investigate spatial and temporal correlations of weather conditions. The spatial resolution is 55 km with 17 varying thickness (exponentially increasing) of horizontal layers. The model temporal resolution is 10 min. The tomography retrieval is based on AWATOS2 model (Perler et al., 2011) that is fed with observations from all GPSNet stations; no external data were introduced. The AWATOS2 uses vertical and horizontal constraints as well as 3D spline parametrisation of troposphere refractivity. The results of tomography processing, compared against radiosonde profiles, agree within 8 ppm (mm/km) in terms of wet refractivity (Manning et al., 2012). 3.3. ZHD and ZWD calculations In order to compute the ZTD, meteorological parameters were converted to the troposphere impacts on GNSS signal using Saastamoinen model (Saastamoinen et al., 1972). The MSL heights from AWS station measurements are converted to the elipsoidal heights by calculating the geoidal undulation N obtained from a global gravimetric geoid model (Boehm et al., 2007). Pressure values are then corrected for the difference Δh Table 1 The validation results for three severe weather events (2010STORM —severe storms [864 reference CODE observations, 144 reference RS observations], 2011FLOOD —flooding [720 reference CODE observations, 192 reference RS observations], 2012FLASH—flash flooding [720 reference CODE observations, 192 reference RS observations]) in terms of ZTD and IWV (both in [mm]). Three data processing strategies were investigated: SHORT —shortest baseline network solution, LNG —longest baseline network solution, and PPP —ZD solution. The missing data for 2011FLOODLNG and 2012FLASHLNG PRED is a consequence of limited number of common observations on two sides of long baselines. PREC term refers to the precise orbits and clocks, whereas the PRED points out to the predicted orbits and real-time clocks. CODE REF [mm] PREC RS REF [mm] PRED PREC PRED ZTD 2010STORMSHORT 2011FLOODSHORT 2012FLASHSHORT 2010STORMLNG 2011FLOODLNG 2012FLASHLNG 2010SSTORMPPP 2011FLOODPPP 2012FLASHPPP μZTD −0.6 −2.3 0 −3.3 – −7.1 −1.7 −1.1 −0.6 σZTD 9.0 7.4 7.6 7.8 – 7.4 10.9 7.2 6.8 μZTD −1.1 −2.4 0 0.8 – – −3.5 −2.7 −0.5 σZTD 8.5 7.7 7.6 13.3 – – 15.4 16.4 13.5 μZTD −15.2 −14.2 −10.1 −12.2 – −21.0 −10.2 −6.4 −4.3 μZTD −9.1 −13.1 −9.8 −12.2 – – −7.2 −9.2 −12.4 IWV 2010STORMSHORT 2011FLOODSHORT 2012FLASHSHORT 2010STORMPPP 2011FLOODPPP 2012FLASHPPP μIWV −2.5 −2.2 −2.2 −1.8 −2.1 −2.0 σIWV 3.1 2.8 2.6 2.0 2.1 1.8 μIWV −1.9 −2.2 −2.5 −1.2 −1.1 −2.4 σIWV 2.6 2.6 2.7 3.5 2.5 10 μIWV −2.7 −2.4 −2.4 −1.7 −0.8 −1.3 μIWV −1.9 −2.2 −2.4 −1.3 −1.1 −2.5 422 between AWS sensor elevation and GNSS antenna elevation using the following equation (Berg, 1948): 5:225 pAWS@GNSS ¼ pAWS
ð1−0:0000226
ðΔhÞÞ : ð19Þ 40 30 20 STAR PKVL MOBS KYNE GEEL EPSO CRES CLAC 10 BALL Conversion between ZTD and IWV requires two types of information: 1) pressure values to remove the dry delay from δZTD normal conditions δZTD storm conditions 50 BACC 3.4. IWV calculation 60 APOL The WV and temperature are not interpolated to the location of the receiver. The assumption made in this study is that 50 km radius weather conditions do not change. The radiosonde pressure observations pRAOBS (2D data) are processed (Eq. (19)) using a procedure similar to the point observations pAWS. Therefore, the ground pressure reading from the radiosonde is interpolated to the height of the receiver. The pressure value is an input for ZHD calculation using the Saastamoinen model (Saastamoinen, 1973). Another approach was applied for WV and temperature. These parameters were interpolated to the horizontal layers of fixed height above the ground. There are two assumptions: 1) meteorological conditions (WV and temperature) at the interpolation layer are similar for GNSS signal and for radiosonde, and 2) the troposphere above the highest layer (in this study 16 km) has negligible impacts on GNSS signal propagation. In order to obtain ZWD, the wet refractivity Nv (Eq. (2)) has been calculated at each horizontal layer and multiplied by the separation distance between layers. The NWP model pressure outputs (3D data) are interpolated to the station locations using the method developed by Bosy et al. (2010) that takes interpolation values from the eight closest NWP model nodes (4 above and 4 below the GNSS station). The pressure value is an input for the ZHD calculation using the Saastamoinen model (Saastamoinen, 1973). The temperature and WV values are interpolated to the 3D grid, similar to the tomography model structure (Rohm and Bosy, 2011). The model horizontal resolution is 100 km, while its vertical resolution is 300 m in the lower troposphere to 2 km near the tropopouse. In each column above the receiver, ZWD is calculated in a similar way to that employed for radiosonde data. The coarse 3D grid horizontal resolution was chosen to account for radiosonde GNSS receiver distance and satisfy the assumption of similar weather conditions on both ends of this baseline. δZTD / δT [mm/h] W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 Fig. 4. The troposphere rate of change in normal and severe weather conditions for CORS stations affected by the storm. the ZTD to retrieve the ZWD (Eq. (6)) and 2) factor Q (Eq. (7)) that translates ZWD (a signal dependent value) to IWV(a water vapour dependent value). The ZWD value is obtained using the ZHD derived from the NWP model outputs according to the methodology presented in the subsection above. To derive factor Q, the water vapour temperature gradient is required (Kleijer, 2004). Usually a global value for mean temperature is employed (Schüler, 2001). It is, however, advised to use local troposphere parameters if available. In this study NWP derived Q factor has been applied, the value of TM is calculated as a quadratic mean of vertical layers (Eq. (8)). 4. Severe weather case studies The severe weather events are often related to excessive amounts of WV present in the troposphere. As such it is anticipated that the WV should be detectable by the GNSS CORS network. Certainly similar studies conducted all around the world confirm that GNSS is capable of detecting passages of all kinds of severe events, such as multi-cell storms (Choy et al., 2013), tropical storms (Seco et al., 2009) and thunderstorms (Seko et al., 2004). Within the scope of this paper three severe weather instances were selected for future assimilation into the ACCESS-A model. Hence, each severe weather event is only shortly introduced as the aim of this study is to provide an optimal troposphere estimates for NWP models. The CORS network used in each case study has a slightly different network extent and the number of receivers. However, there are 3 CORS stations processed within each case study that; 1) are part of troposphere global solution at CODE, 2) are collocated with a radiosonde and 3) are in the NWP model domain. This manuscript only contains figures for the station MOBS, however Australia Fig. 3. The radar image showing the storm passage over Melbourne, 06.03.2010 3:00 UTC to 06.03.2010 4:30 UTC (BoM). 423 W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 ZTD (m m) Fig. 5. The ZTD evolution of weather conditions from GNSS observations, using 127 CORS stations 5.03.2010 to 6.03.2010. numerical statistics are also available for other two stations (HOB2 and TOW2). The overall statistics for each case study are given in Table 1 at the end of this section. The case studies are composed of: 2010 severe storm that hit the Melbourne metropolitan area, 2011 flooding in Victoria and 2012 flash flooding in Victoria, respectively. The first case study will be discussed in great detail and the last two will be used as a supplementary example of GNSS ground-based troposphere estimation robustness. The 2010 March severe storm in Melbourne (BoM, 2010) started on 6 March, initiated by a developing low pressure system and low pressure trough to the west of the south east Australia, developing showers and thunderstorms in the west of the Victoria during the morning. The thunder storms spread eastward and intensified during the first couple of hours of daytime. This thunderstorm developed to the northwest of the Melbourne Metropolitan area (Fig. 3) and moved through the most suburbs of the city from early GPT M. Fig. 6. Comparison of ZTD time series from a number of data sources (GPT, UNB3m, DD GNSS, PPP, AWS, NWP and radiosondes). In GNSS processing precise orbits and clocks were used. 424 W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 PPP Fig. 7. The time series of IWV as derived from number of data sources and GNSS with real-time orbits and clocks. afternoon progressing to the east and then into the southeastern part of Victoria later in the afternoon. The microwave signal propagation conditions were assessed by means of autocorrelation function (20) of wet refractivities Nv: Rðt0; t1Þ ¼ E½ðNvt0 −μ t0 ÞðNvt1 −μ t1 ފ : σ t0
σ t1 ð20Þ The E is an expectation operator, Nvt0 and Nvt1 are the wet refractivity in epochs t0 and t1 respectively, and μt0 and μt1 are the mean values of wet refractivity in the selected epochs (0 and 1). The R(t0, t1) in storm and non-storm conditions was calculated at significant levels of troposphere: boundary layer 0.3 km, mid-troposphere 3.0 km and top troposphere 8.5 km. The decrease of correlation to the threshold value of 0.5 (50%) for severe weather case has been achieved in about 30 min, whereas the normal conditions result in the same threshold value in about 1 h and 15 min. Consequently the δρδt parameter (Fig. 4) characterising the rate of change of the troposphere delay in the Melbourne metropolitan region is also significantly different from normal conditions. The storm and non-storm ZTD observations were separated based on the radar images and the δt was set to 1 h. The relative constraining scheme σrelZTD (13) should adapt to the weather conditions and take the 2 ⋅ σrel ZTD value in case of severe weather conditions (i.e. 40 mm instead of 20 mm). The formal errors of ZTD obtained using standard relative constraints (the same for severe weather and non-severe weather) show an increase of 0.3 − 0.5 mm across all storm-affected stations in all 3 case studies. The GNSS observations during this severe weather event (3rd of March to 9th of March 2010) were collected at all 127 Australian CORS stations. That is 66 stations operating within GPSNet in Victoria, 41 stations from CORSnet-NSW and 10 stations operated by Geoscience Australia. The data were processed according to the strategy described in Section 2. The last processing step covers normal equation inversion and R estimation of unknowns, in Bernese GPS Software ADDNEQ2 module (Dach et al., 2007) is used to manipulate and invert set of normal equations. ZTD parameters were estimated with an hourly time resolution in the final ADDNEQ2 run using minimum constraints imposed over the translation parameters of CORS stations' weekly coordinates and troposphere relative constraints. The severe meteorological conditions (BoM, 2010) resulted in rapid changes of the ZTD parameter. Time derivative of ZTD evolution is shown in Fig. 5. The figure presents the differences between ZTD estimates at starting time 10:00 UTC on the 5th of March 2010 and following epochs (every 3 h) up to 10:00 UTC on the 6th of March 2010. The cold front approaching from the west is sweeping the moist and warm air in the advancing edge of the front, creating highly unstable conditions and developing strong thunderstorms. The point GNSS ZTD observations show (Fig. 6) signature of highly unusual troposphere conditions. The typical point observations from station MOBS show a drastic increase of ZTD while moist air was approaching from the north and then a sudden drop of ZTD just before the storm (Fig. 5 second row, first and second from the left), followed by high and variable ZTD values. To assess the amount of WV that has played an important role in the formation of this severe weather event ZTD observations have been converted to the IWV point observations (Fig. 7) by applying formulas (6) and (7). The Q parameter has been derived from NWP output profiles. By comparing Fig. 6 with Fig. 7 it is clear that both time series closely match. This implies that this is a purely WV induced phenomenon, without large pressure variations. Detailed studies of ZTD/IWV response to severe weather (Champollion et al., 2004; Brenot et al., 2006; Choy et al., 2013), including this one, show that before catastrophic rainfalls IWV increases to the unnoted values of 40–50 mm (Fig. 7 64 DOY). As the severe weather approaches/develops, high variability of retrieved parameters is observed (Fig. 7 DOY 65–66), followed by drop in the total value and significant IWV oscillation as the system dissipates (Fig. 7 DOY 67). W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 The final part of this study is to assess the impact of real time clocks and orbits on troposphere estimates. The effect is visible especially in the PPP data, which is a direct effect of the zero-difference observations solutions, and therefore no reduction of satellite clocks and orbit errors is possible. The impact of real-time satellite clocks and orbits on DD processing is negligible. Table 1 shows mean accuracy measures as indicated by comparing with CODE reference (CODE REF) data and radiosonde measurements (RS REF). The obtained solutions in terms of ZTD/IWV quality is in line with works of other authors (e.g. Karabatić et al. (2011), Jade et al. (2005), Gutman and Benjamin, (2001), and Tregoning et al. (1998)), the expected ZTD retrieval quality should be on the level of 10–15 mm. The biases and standard deviations are in cm-level when compared with radiosondes data or reference ZTD estimates. It is evident that in all case studies the SHORT baseline DD network solution performs best, and the decrease of accuracy between the precise (PREC) and the predicted (PRED) IGS orbits is actually negligible. However this impact starts to be significant in the case of PPP observations, both bias and standard deviation increases. The estimates start to be noisier. The IWV inherits this property which leads to a decreased accuracy between the precise products and real time products. Surprisingly the bias of IWV is larger than ZTD, we believe that this is attributed to the pressure, temperature and water vapour information from NWP model that are introduced to retrieve the IWV. The data assimilation process should not be affected by this bias since it is relative to the GNSS CODE data and not to the radiosonde data. Overall accuracy, regardless of the data processing strategy, is well within the required accuracy of 5 mm with respect to the radiosondes. 5. Summary and conclusion This study is an initial part of a large project aimed at utilising GNSS ground-based troposphere estimation for enhancing NWP models forecasts in Australia. Other applications such as nowcasting using GNSS ZTD data in Australia were discussed in (Choy et al., 2013). It is anticipated that the ground based GNSS data will contribute to increased to 80% (from 40%) improvement in the intensive rainfall forecast (de Haan, 2013). This paper is focused on the GNSS ZTD and IWV processing in challenging weather conditions that caused the increase of ZTD estimation error and signal noise. Several GNSS signal processing strategies have been investigated including; 1) DD (short baselines network solution, long baselines network solution, baseline-by-baseline solution —run but not considered) and 2) PPP (float ambiguity). The final ZTD/IWV product for data assimilation in NWP is going to be delivered to the BoM in near-real time. Therefore the impact of the predicted orbits and clocks on the estimated ZTD values needs to be quantified. Broadly speaking using DD or PPP in conjunction with precise IGS orbits and clocks is of an acceptable quality (IWV error below 3 mm, ZTD error below 15 mm), even in challenging weather conditions. It has to be noted however, that the DD estimates with short baseline formulation are slightly more consistent with reference values than PPP and other strategies considered. The longest baseline strategy that was considered to be the most appropriate one (decorellation of height and 425 troposphere conditions) for lower atmosphere, was not performing exceptionally well. The standard short baseline strategy with network solution was proved to be the most reliable one in the DD strategies family. The PPP approach has confirmed to be robust and reliable with high precision clocks and orbits. However, PPP solution with real-time clocks and predicted orbits is not performing as accurate as the DD solution. From our experiments we concluded that (at least for investigated case studies): 1) severe weather conditions do not interfere with the ability of the CORS network to provide high quality troposphere estimates for NWP data assimilation, 2) the DD and PPP strategies are of equal accuracy in terms of ZTD estimation while using precise orbits and clocks, and 3) in the near real time, and the real time mode the DD solution will perform better. Acknowledgements This research was supported by the Australian Space Research Program project endorsed to a research consortium led by K. Zhang at RMIT University. The GPS ground-based data was retrieved from the GPSNet (owned and operated by the Department of Environment and Primary Industries, Victoria), CORSnet-NSW (owned and operated by the Department of Land and Property Information, New South Wales), Asia-Pacific Reference Frame (operated by Geoscience Australia) and SmartNET Australia. The meteorological data were kindly provided by the Australian Bureau of Meteorology. Fruitful discussions with B. Carter and J. Bosy are kindly acknowledged. References Barlag, S., de Haan, S., Offiler, D., 2004. TOUGH Targeting Optimal Use of GPS Humidity data in meteorology. Technical Report. KNMI, MetOffice. Berg, H., 1948. Allgemeine Meteorologie. Dümmler's Verlag, Bonn. Boehm, J., Niell, A., Tregoning, P., Schuh, H., 2006. Global Mapping Function (GMF): a new empirical mapping function based on numerical weather model data. Geophys. Res. Lett. 15–26. http://dx.doi.org/ 10.1029/2005GL02554. Boehm, J., Heinkelmann, R., Schuh, H., 2007. Short note: a global model of pressure and temperature for geodetic applications. J. Geod. 81, 679–683. http://dx.doi.org/10.1007/s00190-007-0135-3. BoM, 2010. Monthly Weather Review. Victoria. March 2010. Technical Report. Australian Government.Bureau of Meteorology. Bosy, J., Figurski, M., Wielgosz, P., 2003. A strategy for GPS data processing in a precise local network during high solar activity. GPS Solutions 7, 120–129. Bosy, J., Kontny, B., Borkowski, A., 2009. IGS/EPN reference frame realization in local GPS networks. In: Drewes, H. (Ed.), Geodetic Reference Frames. International Association of Geodesy Symposia, vol. 134. Springer, Berlin Heidelberg, pp. 197–203. Bosy, J., Rohm, W., Borkowski, A., Figurski, M., Kroszczyñski, K., 2010. Integration and verification of meteorological observations and NWP model data for the local GNSS tomography. Atmos. Res. 96, 522–530. http://dx.doi.org/10.1016/j.atmosres.2009.12.012. Brenot, H., Ducrocq, V., Walpersdorf, A., Champollion, C., Caumont, O., 2006. GPS zenith delay sensitivity evaluated from high-resolution numerical weather prediction simulations of the 8–9 September 2002 flash flood over southeastern France. J. Geophys. Res. - Atmos. 111. Champollion, C., Masson, F., Van Baelen, J., Walpersdorf, A., Chéry, J., Doerflinger, E., 2004. GPS monitoring of the tropospheric water vapor distribution and variation during the 9 September 2002 torrential precipitation episode in the Cévennes (southern france). J. Geophys. Res. - Atmos. 109. Choy, S., Wang, C., Zhang, K., Kuleshov, Y., 2013. GPS sensing of precipitable water vapour during the march 2010 Melbourne storm. Adv. Space Res. 52, 1688–1699. http://dx.doi.org/10.1016/j.asr.2013.08.004. Christian, M.H., Wang, C., Angkasa, K., Gritton, K., 2004. Estimation of Microwave Power Margin Losses Due to Earth's Atmosphere and 426 W. Rohm et al. / Atmospheric Research 138 (2014) 414–426 Weather in the Frequency Range of 3–30 GHz. Technical Report.Jet Propulsion Labolatory, California Institute of Technology. Collins, P., Bisnath, S., Lahaye, F., Héroux, P., 2010. Undifferenced GPS ambiguity resolution using the decoupled clock model and ambiguity datum fixing. Navig. 57, 123. Dach, R., Hugentobler, U., Fridez, P., Meindl, M., 2007. Bernese GPS Software Version 5.0. Astronomical Institute, University of Bern. de Haan, S., 2013. Assimilation of GNSS ZTD and radar radial velocity for the benefit of very-short-range regional weather forecasts. Q. J. R. Meteorol. Soc. http://dx.doi.org/10.1002/qj.2087 (n/a–n/a). Douša, J., 2010. The impact of errors in predicted GPS orbits on zenith troposphere delay estimation. GPS Solutions 14, 229–239. http:// dx.doi.org/10.1007/s10291-009-0138-z. Gutman, S.I., Benjamin, S.G., 2001. The role of ground-based GPS meteorological observations in numerical weather prediction. GPS Solutions 4, 16–24. http://dx.doi.org/10.1007/PL00012860. Herring, T.A., King, R.W., McClusky, S.C., 2010. GAMIT Reference Manual Reference Manual GPS Analysis at MIT. Department of Earth, Atmospheric, and Planetary Sciences.Institute of Technology, Massachusetts. Hofmann-Wellenhof, B., Lichtenegger, H., Wasle, E., 2008. GNSS—Global Navigation Satellite Systems: GPS, GLONASS, Galileo, and More. Springer. Jade, S., Vijayan, M., Gaur, V., Prabhu, T.P., Sahu, S., 2005. Estimates of precipitable water vapour from GPS data over the Indian subcontinent. J. Atmos. Sol. Terr. Phys. 67, 623–635. http://dx.doi.org/10.1016/j.jastp.2004.12.010. Karabatić, A., Weber, R., Haiden, T., 2011. Near real-time estimation of tropospheric water vapour content from ground based GNSS data and its potential contribution to weather now-casting in Austria. Adv. Space Res. 47, 1691–1703. Kleijer, F., 2004. Troposphere Modeling and Filtering for Precise GPS Leveling. Ph.D. thesis. Technische Universiteit Delft. Kouba, J., 2009. A guide to using International GNSS Service (IGS) products. Unpublished (available electronically at http://igscb.jpl.nasa.gov/igscb/ resource/pubs/UsingIGSProductsVer21.pdf). Kubo, Y., Tanaka, H., Ohashi, M., Sugimoto, S., 2012. Long baseline GNSS relative positioning with estimating ionospheric and tropospheric delays and their gradients. Int. J. Innov. Comput. I 8, 2375–2388. Leandro, R., Langley, R., Santos, M., 2008. UNB3m_pack: a neutral atmosphere delay package for GNSS. GPS Solutions 12, 65–70. Li, X., Ge, M., DouÅ¡a, J., Wickert, J., 2013. Real-time precise point positioning regional augmentation for large GPS reference networks. GPS Solutions 1–11. http://dx.doi.org/10.1007/s10291-013-0310-3. Manning, T., Zhang, K., Rohm, W., Choy, S., Hurter, F., 2012. Detecting severe weather using GPS tomography: an Australian case study. J. Global Pos. Syst. 11, 58–70. http://dx.doi.org/10.5081/jgps.11.1.58. Marshall, J., Xiao, Y., Norman, R., Zhang, K., Rea, A., Cucurull, L., Seecamp, R., Steinle, P., Puri, K., Le, T., 2010. The beneficial impact of radio occultation observations on Australian region forecasts. Aust. Meteorol. Ocean J. 60, 121–125. Mendes, V., Langley, R., 1999. Tropospheric zenith delay prediction accuracy for high precision GPS positioning and navigation. Navig. 46, 25–34. Misra, P., Enge, P., 2001. Global Positioning System: Signals, Measurements, and Performance. Ganga-Jamuna Press. Owens, J., 1967. Optical refractive index of air: dependence on pressure, temperature and composition. Appl. Opt. 6, 51–59. Perler, D., Geiger, A., Hurter, F., 2011. 4D GPS water vapor tomography: new parameterized approaches. J. Geod. 85, 539–550. http://dx.doi.org/ 10.1007/s00190-011-0454-2. Rohm, W., Bosy, J., 2011. The verification of GNSS tropospheric tomography model in a mountainous area. Adv. Space Res. 47, 1721–1730. http:// dx.doi.org/10.1016/j.asr.2010.04.017. Saastamoinen, J., 1973. Contribution to the theory of atmospheric refraction, in three parts. B. Géod. 105 (106), 107 (279–298, 383–397, 13–34). Saastamoinen, J., et al., 1972. Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites. Geophys. Monogr. Ser. 15, 247–251. Schüler, T., 2001. On Ground-based GPS Tropospheric Delay Estimation. Ph.D. thesis Universitat der Bundeswehr Munchen. Schüler, T., 2006. Impact of systematic errors on precise long-baseline kinematic GPS positioning. GPS Solutions 10, 108–125. Seco, A., González, P., Ramírez, F., García, R., Prieto, E., Yagüe, C., Fernández, J., 2009. GPS monitoring of the tropical storm delta along the Canary Islands Track, November 28–29, 2005. Pure Appl. Geophys. 166, 1519–1531. Seko, H., Nakamura, H., Shoji, Y., Iwabuchi, T., 2004. The meso-γ scale Water vapor distribution associated with a thunderstorm calculated from a dense networks of GPS receivers (3. application of GPS data to atmospheric science). J. Meteorol. Soc. Jpn 82, 569–586. Solheim, F., Vivekanandan, J., Ware, R., Rocken, C., 1999. Propagation delays induced in GPS signals by dry air, water vapor, hydrometeors, and other particulates. J. Geophys. Res. 104, 9663–9670. Strang, G., Borre, K., 1997. Linear Algebra, Geodesy, and GPS. Wellesley Cambridge Pr. Thayer, G.D., 1974. An improved equation for the radio refractive index of air. Radio Sci. 9, 803–807. Tregoning, P., Boers, R., O'Brien, D., Hendy, M., 1998. Accuracy of absolute precipitable water vapor estimates from GPS observations. J. Geophys. Res. - Atmos. 103, 28701–28710. http://dx.doi.org/10.1029/98JD02516.