Review of the state of the art and future prospects of the ground-based GNSS meteorology in Europe

Global navigation satellite systems (GNSSs) have revolutionised positioning, navigation, and timing, becoming a common part of our everyday life. Aside from these well-known civilian and commercial applications, GNSS is now an established atmospheric observing system, which can accurately sense water vapour, the most abundant greenhouse gas, accounting for 60–70 % of atmospheric warming. In Europe, the application of GNSS in meteorology started roughly two decades ago, and today it is a well-established field in both research and operation. This review covers the state of the art in GNSS meteorology in Europe. The advances in GNSS processing for derivation of tropospheric products, application of GNSS tropospheric products in operational weather prediction and application of GNSS tropospheric products for climate monitoring are discussed. The GNSS processing techniques and tropospheric products are reviewed. A summary of the use of the products for validation and impact studies with operational numerical weather prediction (NWP) models as well as very short weather prediction (nowcasting) case studies is given. Climate research with GNSSs is an emerging field of research, but the studies so far have been limited to comparison with climate models and derivation of trends. More than 15 years of GNSS meteorology in Europe has already achieved outstanding cooperation between the atmospheric and geodetic communities. It is now feasible to develop next-generation GNSS tropospheric products and applications that can enhance the quality of weather forecasts and climate monitoring. This work is carried out within COST Action ES1206 advanced global navigation satellite systems tropospheric products for monitoring severe weather events and climate (GNSS4SWEC, http://gnss4swec.knmi. nl).


Introduction
Atmospheric water vapour has a complex life cycle in the troposphere, including vertical and horizontal transport, mixing, condensation, precipitation, and evaporation.Due to its high temporal variation (more than 50 % within a few hours, Johansson et al., 1998) and its complex distribution, linked to its relationship with atmospheric dynamics and the role Published by Copernicus Publications on behalf of the European Geosciences Union. of phase changes, the understanding of tropospheric water vapour is essential to climate science and numerical weather prediction (NWP).However, while being a critical atmospheric observation parameter, it is also very demanding to observe.Since the 1950s, the standard technique for measuring atmospheric water vapour has been the radiosonde.
While radiosondes provide invaluable atmospheric observations, they do have some inherent limitations, primarily limited temporal and spatial resolution (soundings are typically performed once or twice a day and the spatial coverage of the network in Europe is 250 km or greater).Tralli and Lichten (1990) first proposed using the global positioning system (GPS) for atmospheric sounding.As GPS signals travel through the atmosphere, their propagation is affected by atmospheric constituents, in particular water vapour.The technique was first called GPS meteorology and later renamed to global navigation satellite system (GNSS) meteorology.GNSS refers to any satellite constellation with global coverage used for positioning, navigation, and timing.Currently, two constellations are in full operation: GPS (USA, NAVSTAR) and GLONASS (Russia).Two other systems are currently under development: Galileo (ESA/EC, with 12 satellites in orbit out of an operational constellation of 30) and BeiDou (China, with 17 satellites in orbit out of an operational constellation of 35).In this paper we will use the name "GPS meteorology" when dealing with the stand-alone GPS and the name "GNSS meteorology" when referring to a system other than this GNSS or in a more general context.The establishment of GPS meteorology as an operational atmospheric sounding technique in Europe was the focus of a number of major European projects (Fig. 1).Collaborative activities aiming at the exploitation of GPS observations for monitoring atmospheric water vapour started with the European Commission (EC) 4th Framework Program (FP) projects WAVEFRONT (GPS WAter Vapour Experiment For Regional Operational Network Trials) and MAGIC (Meteorological Applications of GPS Integrated Column Water Vapour Measurements in the western Mediterranean).MAGIC combined the efforts of GPS network operators in France, Italy, and Spain to build a GPS demonstration network for weather forecasting in the western Mediterranean area (Haase et al., 2001).In the follow-up project, EC COST Action 716 (COST-716), 15 European countries participated with the overall aim of providing a European demonstration campaign for GPS meteorology.Operational provision of ground-based GPS zenith total delay (ZTD) products started in 2001 (Elgered et al., 2005).In 2003, MAGIC and COST-716 were followed by the EC 5th FP scientific project TOUGH (Targeting Optimal Use of GPS Humidity Measurements in Meteorology), and in April 2005, the EIG EUMETNET GPS Water Vapour Programme (E-GVAP, http://egvap.dmi.dk) was established to transform GPS meteorology from an R&D activity to a fully operational service across Europe.Currently, 18 National Meteorological and Hydrological Services (NMHSs) collaborate in E-GVAP with 17 GPS analysis centres (ACs), which collect and process GPS data from over 1800 European ground-based stations (Fig. 2).
The present state of the art of GNSS meteorology would have been impossible without the products provided by the Global Geodetic Observing System (GGOS, http://www.ggos.org).GGOS products are provided by the services of the International Association of Geodesy (IAG, http:// www.iag-aig.org).Of primary importance is the International GNSS Service (IGS, http://www.igs.org).At present, the IGS provides a wide range of satellite (orbit and clocks) and tropospheric products essential for GNSS meteorology.
The aim of this paper is to assess the current state of the art of GNSS meteorology in Europe.In Sect.2, the development of GNSS tropospheric product processing techniques are reviewed.Section 3 summarises the use of GNSS tropospheric products for operational weather prediction applications like NWP, nowcasting, and monitoring severe weather.Section 4 presents the application of GNSS tropospheric products in climate monitoring.The outlook and future work are given in Sect. 5. 2 GNSS tropospheric processing 2.1 GNSS observations and the atmosphere GNSS signals are delayed and bent when propagating through the atmosphere.The upper part of the atmosphere is a dispersive medium, the influence of which can be more or less eliminated by combining observations from the two GNSS L-band frequencies in the range from 1.16 to 1.61 GHz.The effect of the lower neutral part of the atmosphere, however, cannot be eliminated in a similar way because it is a non-dispersive medium.The major effect of the neutral atmosphere occurs in the troposphere.In this paper we refer to this effect as tropospheric delay.In accurate GNSS applications, and specifically for our application, the tropospheric delay at each station is estimated, together with other parameters, in the GNSS data analysis.
A common model for the total slant path delay from the GNSS satellite to the receiver on the surface of the Earth (see e.g.Teke et al., 2013) is as follows: ρ sat rec (e, A) = mf h (e) ZHD + mf w (e) ZWD (1) where e and A are the elevation and the azimuth angle towards a specific satellite, ZHD and ZWD are the zenith hydrostatic and zenith wet delays expressed in units of metres, mf h , mf w , mf g are the hydrostatic, wet, and gradient mapping functions, and G N and G E are the components of linear horizontal gradients.The first two terms on the righthand side represent models assuming tropospheric symmetry while the last term may be added in order to estimate a firstorder asymmetry in terms of a linear horizontal gradient.
The tropospheric zenith total delay (ZTD) characterises the effect of the troposphere which delays GNSS signal were they come from the zenith.The mean delay is typically about 2.4 m from a site at mean sea level.The ZTD is the sum of the ZHD and the ZWD or, alternatively, a sum of an a priori and an estimated tropospheric delay correction.A separation of hydrostatic and non-hydrostatic components was proposed by Davis et al. (1985).The latter is often simply referred to as the ZWD although it represents a non-hydrostatic component complementary to the hydrostatic one.The elevation angle dependencies of ZHD and ZWD are modelled by individual mapping functions with coefficients calculated using climatology, meteorologic observations or NWP output.
The definitions of ZHD and the ZWD originates from the integrals describing the delay and bending angle of a signal propagating through the atmosphere with a varying refractivity defined by meteorological parameters.A model for the ZHD was developed by Henriksen et al. (1972): where where k 1 is a refractivity constant, K Pa −1 , R d is the specific gas constant of dry air (287.058J kg −1 K −1 ), g ϕ,H is the gravitational acceleration in metres per second squared at  the latitude ϕ and at the station height H , in metres above the geoid, p is the atmospheric pressure in Pa, at the height H in metres.
Since ZWD is mainly a function of partial water vapour pressure (e) and the temperature (T ), it can be calculated with a numerical integration through a full atmospheric profile using these two meteorological parameters together with the two refractivity constants k 2 [K hPa −1 ] and The above description summaries the relations between tropospheric parameters used in the GNSS data analysis and the meteorological parameters of interest.A synergy between GNSSs and meteorological observations thus offer broad range of applications.Firstly, the capability of estimating tropospheric parameters from GNSS with a high spatial and temporal resolution and low latency is of great interest in meteorology.Secondly, since the long-term GNSS observations are available globally, monitoring trends in the atmospheric water vapour is also feasible.Thirdly, meteorological or climatologic data are useful when trying to optimise the models used in GNSS data processing.Fourthly, meteorological or climatologic data can be useful to support high-rate realtime kinematic applications in order to reduce tropospheric errors or, more generally, to avoid high correlations between estimated height coordinates.Within GNSS tropospheric processing, it is common to distinguish products according to their availability (Table 1): (1) final products (Final) based on daily observations and final precise post-processed products, (2) near-real-time products (NRT) based on hourly observations and ultra-rapid precise products, (3) real-time products (RT) based on real-time observations and real-time precise products.The final products (column 2 in Table 1) have the best possible precision and accuracy, but are available with a delay of between 1 day and 2 weeks after the time of observation.Due to their high quality, the final products are used as reference products.The NRT and RT (columns 3 and 4 in Table 1) product accuracy is compromised due to the latency requirements implied by the application.For example, in order to be used in NWP models, the NRT products which are updated hourly have a requested latency of 90 min after time of observation.The RT products have a typical latency of 5-10 min after time of observation and are updated sub-hourly.
The potential of GPS observations for tropospheric monitoring was suggested by Tralli and Lichten (1990) and Bevis et al. (1992), but it only concerned Final Products.The development of NRT GPS data processing suitable for operational weather prediction was, however, not possible until 1999 mainly due to (1) the lack of an hourly data flow from the permanent GPS stations and (2) the limited quality of predicted precise satellite orbits and clocks.First operational NRT tropospheric products (ZTD) for Europe were demonstrated in the early 2000s (Gendt et al., 2001;Dick et al., 2001;Douša, 2001a, b) during the COST-716 benchmark campaign.After this successful benchmark campaign, a demonstration campaign was initiated in 2001 (van der Marel et al., 2004).The demonstration campaign continued without interruption during the TOUGH and E-GVAP I, II, III, and IV projects (see Fig. 1).The main purpose of the E-GVAP project is to collect and distribute GNSS-derived ZTDs for usage in operational weather forecasting in NRT (within 90 min), to monitor data quality and attempt to expand the GNSS meteorological network in a collaboration between meteorology and geodesy.

Advanced GNSS processing techniques for NRT and RT products
The current state of the art in GNSS data processing for operational meteorology provides hourly updated NRT ZTD estimates, using mainly data from only the GPS satellite constellation.The majority of analysis centres (ACs) utilise a network solution to eliminate GNSS receiver and satellite clock errors in double-difference observations.Two GNSS ACs provide an alternative approach based on the precise point positioning (PPP) technique (Zumberge et al., 1997), which exploits original data without differencing.The main reason for applying the network approach was the lack of ultrafast precise satellite clock products required for PPP (see Table 2).The ACs GFZ (German Research Centre for Geosciences, Gendt et al., 2004)    (2) an autonomous, distributed PPP ZTD processing for each station individually.Significant improvements in GNSS data processing (including models and products), as well as enhanced GNSS infrastructure, lead to new challenges.The real-time data and product flows have been developed and standardised, and many GNSS stations from around the world now provide real-time data streams, which are used for generating realtime precise products such as precise satellite clock corrections.In 2013, the International GNSS Service (IGS, Dow et al., 2009) launched its real-time service (Caissy et al., 2012), providing precise satellite orbit and clock corrections in support of ultra-fast or real-time PPP processing (Douša and Vaclavovic, 2014;Li et al., 2014a;Yuan et al., 2014).Tropospheric products derived from real-time PPP processing are of high interest for NWP models and for nowcasting of severe weather events.
Along with the standard estimation of ZTD, several ACs monitor the tropospheric asymmetry by estimating tropospheric horizontal gradients.In addition, retrieval of the slant path delay (SPD), i.e. the total tropospheric delay in the direction from the receiver to a visible satellite, can provide additional information about the atmospheric structure (see Fig. 3).One important application of SPDs is their use in tomographic reconstruction of 3-D water vapour fields (see Sect. 3.3).However, tomographic reconstruction depends to a large extent on the quality of the SPD estimates.The PPP processing method for producing SPDs is preferred due to its efficiency, namely (1) easy support of additional estimated parameters in the distributed processing system, and (2) the capability of directly providing undifferenced slant delays from the PPP.However, this strategy requires precise satellite orbits, clocks products, and the most accurate models.
The processing of individual SPDs is currently an active field of research and different approaches are considered.One strategy is to rely on the combination of the wellestablished production of ZTDs and horizontal gradients to map this information to the direction of the GNSS signal path.Several variations are used: (1) mapping of the zenith delay on the slant axis (Notarpietro et al., 2011), (2) mapping of the zenith delay with gradients (Brenot et al., 2014;Kačmařík and Rapant, 2012), in some cases separating the dry and wet parts and using the appropriate mapping function (Champollion et al., 2009) or (3) mapping of the zenith delay with gradients and adding post-fit residuals (Bender et al., 2011b).Another aspect is the separation of the SPD to slant hydrostatic delay (SHD) and slant wet delay (SWD).The GNSS processing cannot distinguish between the two parts and a combination of meteorological observations and empirical models is required to estimate the ZWD and the corresponding SWD (Bevis et al., 1994).This can either be done within the SPD processing or as an independent step after the GNSS processing.Considerable errors can be introduced within the separation step as the SWD is only a small part of the SPD (0-10 %) and small errors in the SHD can lead to severe errors in the SWD, especially for dry periods.All of these approaches were used to process SPDs for GNSS tomography experiments, and realistic humidity fields could be obtained with any of the SPD data sets.Even in direct comparison, it is very difficult to identify the optimum strategy (Kačmařík et al., 2012).
Emerging new satellite systems and their signals might be used for strengthening ZTD solutions and for improving the estimates of atmospheric asymmetry in the vicinity of the GNSS receiver.Although the GLONASS system is still not officially in operation, data are already being collected and used.Since 2008, the IGS has provided ultra-rapid precise orbits for GLONASS satellites with contributions from at least four European ACs with products being sufficiently robust for generating operational multi-GNSS (GPS and GLONASS) NRT ZTD solutions (Douša, 2010).The initial inconsistency between stand-alone GPS and GLONASS ZTD solutions disappeared when the IGS08 antenna phase centre models were adopted (Dach et al., 2011b).Li et al. (2014b) reports on the first results of tropospheric monitoring with the Chinese BeiDou GNSS system.Initial problems still need to be resolved, namely enhanced processing software including strategies for optimal multi-signal and multifrequency processing, support of precise orbit and clock products for the BeiDou satellites, and improved models for receiver and satellite phase antenna offsets and variations.
In 2010, the first operational global NRT ZTD estimates were provided by ACs at the UK Met Office and at GOP (Geodetic Observatory Pecny) with data from the GOP solution having been assimilated into the Met Office global NWP model since 2012 (Douša and Bennitt, 2013).Both systems apply the standard double-difference ZTD processing scheme.However, in the future, the PPP technique can be applied in order to support an unlimited number of stations, with an almost homogeneous ZTD quality expected with some deviations due to regional quality variations in the global precise products.
The use of meteorological data as input in GNSS processing has been limited.Climatological data were usually preferred for the utilisation in GNSS analyses mainly for (1) developing mapping functions (which represent the tropospheric path delays at any elevation derived from the tropospheric delay at zenith above the GNSS antenna) and ( 2) separating hydrostatic and wet tropospheric components of the ZTD.In the last decade, the increasing demands of accuracy for positioning and navigation, increasing resolution of NWP models, and the direct use of actual meteorological data in space-based geodetic techniques has been explored.For example, the Vienna Mapping Function (VMF) concept provides a support for mapping functions derived from NWP models by applying the ray-tracing technique (Böhm and Schuh, 2004;Urquhart et al., 2014).The impact of atmospheric loading based on actual meteorological data has also been evaluated by various ACs (e.g.Dach et al., 2011a).Future services supporting GNSS with meteorological data as input will further improve the temporal and spatial data resolution and will include modelling of atmospheric asymmetry (e.g.Böhm and Schuh, 2007).Other emerging fields requiring exploitation of actual meteorological data are real-time precise positioning, navigation, and timing applications.

Reprocessed GNSS tropospheric products
Within EUREF (http://www.euref.eu),tropospheric parameter estimation started in 2001 as a Special Project.In 2008, the Special Project was declared successful and the provision of the tropospheric product became a routine operation with daily delivery of a combined tropospheric solution.The EUREF Permanent Network (EPN, http://epncb.oma.be) is a key geodetic infrastructure coordinating a network of about 250 European GNSS stations.For each station, EUREF provides daily tropospheric estimates at 1-h sampling rate for the period 1996-present.More than 100 EPN stations are co-located with radiosondes, while only a small subset of stations are collocated with very long baseline interferometry (VLBI) radio telescopes.At these collocated sites, routine comparisons are carried out, which can be used for intertechnique evaluation (see Sect. 2.4).Taking benefit of more than 17 years of continuous GNSS observations in Europe, the second reprocessing activity (EPN-repro2) is ongoing using the latest available models and analysis strategy in order to provide an homogeneous GNSS time series suitable for monitoring climate.In 2007, a Memorandum of Understanding between EUREF and EIG EUMETNET was signed to facilitate data exchange and to promote increasing cooperation between the geodetic and meteorological communities in Europe (Pottiaux et al., 2009).
Within the IGS, global GNSS tropospheric products have been produced since 1997.Until 2003, ZTD estimates consisted of a combined product derived from the contribution of the routinely produced ZTDs from several IGS ACs (Gendt, 1998).Due to large inconsistencies among analysis centres and changes over time in the estimation processes, the legacy combination product was abandoned in 2003 and was replaced with a new ZTD product derived from PPP processing by one AC (Byun and Bar-Sever, 2009).In 2008, the IGS initiated the first reprocessing campaign (IGS-repro1, http://acc.igs.org) in order to provide a homogeneous time series of the IGS core products and relevant derived products including tropospheric estimates.The reprocessing was done in three sequential steps: (1) nine ACs delivered global precise products, (2) those contributions were combined to produce the official IGS-repro1 orbit and clock products, and (3) homogeneous tropospheric parameters were calculated using PPP for the period 1995 to 2007 (Bar-Sever, 2012).Using the same processing procedure, the IGS ZTD product was updated for the period 2007 to 2010 and continued operationally until mid-2011, giving access to a homogeneous ZTD series of more than 16 years.Since April 2011, the IGS tropospheric parameters are produced by the United States Naval Observatory (USNO), using PPP processing, but with Bernese software (Byram and Hackman, 2012).In 2013, a second IGS reprocessing campaign (IGS-repro2) was initiated, from which a new reprocessed tropospheric product is expected.
In the last 5 years, reprocessing activities have been ongoing in many national agencies and European ACs (regional and local reprocessing activities).It is to be noted that all reprocessing strategies require well-documented, long-term metadata.Thus, an effort will be required to archive data and metadata from the regional and national densification networks for future reprocessing activities.

Intertechnique comparisons
Intercomparisons are made for a number of reasons.When a new technique becomes available, as in the case of GNSS, it is naturally compared to existing techniques to assess its strengths and weaknesses.Atmospheric water vapour is difficult to measure with high accuracy, so when the GNSS technique appeared in the mid-1990s it was almost immediately also used as a reference, e.g. for evaluating NWP models.
Weather balloons carrying radiosondes vertically up through the atmosphere are the classical meteorological in situ observation technique.In terms of its usage for direct comparison, the main disadvantages of radiosondes are that the measurement direction cannot be controlled (its horizontal path is determined by the wind speed and direction) and that it is not an instantaneous observation (the ascension of the balloon typically takes 1-2 h).Both of these factors must be taken into account when radiosonde data are compared to other remote sensing observations.
Techniques capable of remote sensing atmospheric water vapour can be divided into two categories: (1) differential time of arrival measurements, and (2) emission/absorption measurements.Differential time of arrival techniques were originally developed for positioning, navigation, and timing applications with atmospheric water vapour observation being a spin-off product from space-based geodetic data processing techniques.In addition to GNSSs, there are two other techniques which use signals in the radio band of the electromagnetic spectrum, namely VLBI and Doppler Orbitography Radiopositioning Integrated by Satellite (DORIS).
Microwave radiometry from the ground measures the emission from atmospheric water vapour around the 22 GHz (MWR22) or the 180 GHz (MWR180) emission line.The 22 GHz system is more common since it also works reasonably well in cloudy conditions, being able to correct for the influence of cloud water droplets on the total emission observed by using a second observation frequency.Both systems, however, do not provide any meaningful observation of atmospheric water vapour during rainy conditions.Additional ground-based techniques, which may observe in the visible band, comprise the CIMEL sun photometers from the Aerosol Robotic Network (AERONET), and the Fourier transform infrared spectrometers (FTIRs) observing solar absorption spectra.
The Advanced Microwave Sounding Unit-B (AMSU-B) is a satellite-based microwave instrument that can be used to infer the atmospheric water vapour content from several bands around the 180 GHz emission line from polar orbiting meteorological satellites.In the infrared and visible bands, satellite measurements are made with the Atmospheric Infrared Sounder (AIRS) and the GOME/SCIAMACHY/GOME-2 instruments, respectively.
Table 3 summarises a number of published studies that include the aforementioned instruments and techniques.These, together with references therein, provide an extensive comparison of the results.It is, however, not possible to draw final conclusions regarding the absolute accuracies of the instruments and techniques.All techniques suffer from systematic errors at different timescales meaning that case studies using specific sensors, at specific locations, at different times, will give different comparison results, often presented as biases and root mean square (RMS) differences.For example, an error appearing as a bias in a 2-week long comparison may present itself as a random error over a period of many years.
It is not obvious that any of the techniques listed in Table 3 is superior in terms of accuracy: GNSS, for example, suffers from systematic errors that depend on the elevation cut-off angle used in the processing scheme.By increasing elevation angle, the short-term (random) error increases but the stability over longer time periods is improved as systematic errors are more evident at low elevation angles (Ning and Elgered, 2012).The VLBI technique is not affected by multipath due to the large directivity of the radio telescopes.However, even though VLBI may be the most accurate technique available today, it has a large disadvantage in terms of acquiring operational data as VLBI measurements are only made at tens of sites globally, the observations are not continuous, and each telescope can only measure in one direction at the time.
In conclusion, ground-based GNSS has short-term RMS errors around 3-4 mm in the ZTD, corresponding to 0.4-0.6 kg m −2 in integrated water vapour (IWV).Systematic errors include (1) satellite and receiver instrumentation effects (antenna phase centre offsets and their variations, receiver code, and phase delays), (2) in situ environmental effects (multipath, effects due to snow and ice on the antenna, monumentation), and (3) various modelling approximations such as mapping functions and methods for separation of the ZHD and the ZWD.For NRT applications and weather forecasting, systematic errors are less serious since it is the detection and timing of IWV changes that offers valuable information.For the extreme application of climate monitoring where stability over decades is required, GNSS is still a valuable technique, but here other space geodetic systems such as DORIS, and especially VLBI, can fulfil the role of maintaining stability through collocated observations.

GNSS in NWP and weather forecasting
Comparisons of GNSS tropospheric products to NWP model data, as well as the use of the GNSS data in weather forecasting, have been investigated by many groups.Here we focus mainly on publications which discuss this topic from a European perspective.The review is divided into subsections covering intercomparison of GNSS and NWP model output (Sect.3.1.1),data assimilation (DA), and impact of GNSS data in NWP models (Sect.3.1.2),use of GNSS data by forecasters for nowcasting and severe weather monitor-ing (Sect.3.2), and finally GNSS water vapour tomography (Sect.3.3).

Usage of GNSS tropospheric products in NWP models
While much of the work regarding application of GNSS data in weather forecasting is taking place at the National Meteorological and Hydrological Services (NMHSs), only some of these results are reported in refereed journals.In a few cases, non-reviewed and unpublished research is included in this review in order to accurately present the current state of the art regarding the use of ground-based GNSS data in meteorology.

GNSS tropospheric products comparison with NWP
In Europe a multitude of NWP models are used by the NMHSs.For an overview of the configurations used, see the EUMETNET-SRNWP compilation at http://srnwp.met.hu.The majority of these models are developed in collaboration between NMHSs in model consortia, see e.g.http: //srnwp.met.hu/Consortia/ListofConsortia.html.A few additional, more freely available, NWP models are used at universities and other research institutions.Each NWP model has its own characteristics, resulting in parts of their behaviour being similar from institute to institute.However, NWP models are set up slightly differently to suit local requirements and both the GNSS data processing and NWP models have been improved substantially during the time covered by this review.
For these reasons, biases between NWP and GNSSs change with the model, but they also change with region and with time, similarly to changes in NWP model skills.Most of these variations are gradual, while a few are on/off in nature.Examples include the move to absolute antenna phase centre calibration in GNSS data processing, which changed the GNSS bias relative to all NWP models, and the use of observation operators and height corrections in NWP DA.The observation operators (calculating the GNSS equivalent, for example ZTD, from the NWP model variables) in different NWP models are not identical.In particular, they treat the ZTD contribution above the model top and the corrections for height offsets between GNSS antenna and model orography differently.These differences cause "artificial" systematic ZTD biases between NWP models, adding to the "real" biases resulting from slightly different humidity, pressure, and temperature fields.As long as the biases caused by the observation operator are not (or only weakly) weather dependent, they can effectively be negated by a bias correction prior to DA.However, this is something to be aware of when validating NWP against GNSS ZTD data.Use of GNSS data to validate the quality of NWP diurnal cycle variations in humidity is much less affected by the aforementioned bias changes.It should be noted that much of the work on comparisons cited below is from a time when the main issue was to document that GNSS ZTD estimates were useful to meteorology.A fundamental change took place in 2005, when E-GVAP was established to provide GNSS ZTDs specifically for operational meteorology, and the focus moved toward usage rather than validation.Table 4 summarises comparisons of GNSS tropospheric products and NWP models in Europe.It reports (1) a pronounced seasonal cycle of GPS-NWP model bias and standard deviation (SD) (COSMO and HIRLAM), with the largest errors in summer and smallest in winter; (2) bias and SD decreased with altitude, but bias was almost constant with latitude, while SD decreased with latitude (HIRLAM).In the COSMO NWP model, underestimation of IWV was found to correlate well with significant overestimation of light precipitation, suggesting a coupling between the two components of the hydrological cycle.Also reported is a systematic underestimation of NWP diurnal IWV cycle between 06:00 and 21:00 UTC (COSMO).The largest differences are observed at stations located in mountainous areas and/or near the sea, which reveal differences in model representativeness (ECMWF).Large GNSS-NWP model ZTD biases are found over limited periods of 1-3 days, mostly at coastal stations (ALADIN).

Assimilation of GPS tropospheric products in NWP models -impact studies and long-term experiments
NWP is an initial value problem, the quality of the initial state is of crucial importance to the quality of the forecasts.
The procedure of finding the initial state is called DA.
In most cases, such as in the most commonly used 3-D-Var and 4-D-Var (Var means variational) DA systems, the initial state is obtained by combining observations with an NWP field (a short-term forecast from the previous forecast cycle, the so-called first guess, valid at the time of the data analysis) before the start of the NWP forecast run.This is done in a statistically optimal way.The differences between 3-D-Var and 4-D-Var are that, in the latter, the time evolution within the DA time window (i.e. the period inside in which the observations are selected) is taken into account and that a time sequence of observations from a site can be used, while only one observation is used in 3-D-Var.At the present time, 4-D-Var is mainly used for global models, with wide assimilation time windows (e.g.6-12 h), but in some cases it is also used in regional models with more narrow assimilation windows.
Due to the inclusion of an error covariance matrix for the NWP model and of physical balances in Var assimilation, a given observation will influence several model variables simultaneously and in an extended region around the location of the observation.Each observation is associated with observation errors, including both the instrumental error and the error of representativeness.The latter is caused by the fact that the volume represented by the observation often differs from the resolution of the model.In most cases the observations are more local, i.e. for wind and temperature.Development of a Var-based assimilation system requires manpower and a powerful computer to run in an operational setting.It must be noted that different institutions using the same NWP model often use different approaches for ZTD bias correction as well as for the selection of observations and determination of error statistics, hence the impact will vary from institution to institution.An alternative NWP DA technique is nudging.Instead of first making a separate data analysis, nudging is done by adding to the basic NWP model equations artificial forcing terms, the strength of which depends on the deviation between the NWP and the observations.For the part of the new forecast cycle for which fresh observations are available, one adds or "nudges" in the observations.When the present time has passed (i.e.no more new observations are available), the forecast becomes a truly free forecast.The benefit of nudging is its simplicity, the ease with which a whole time sequence of observations can be used, and the ability to use very recent observations.Some drawbacks of the method are that one can only nudge in observations that correspond to model variables or are somehow transformed to correspond to model variables and the less sophisticated quality control of observations.In practice, nudging GNSS data is done by first converting ZTD to IWV, then nudging humidity up or down in the column in the NWP model, depending on the NWP GNSS offset.
The impact of a new type of observations in NWP can be investigated in many ways.The most stringent is observing system experiments (OSE), in which two simulations are performed, one with and one without (the control run) assimilation of the new observations.It is important to observe which observations were used in the control runs.Often one will only consider the new observations in the development phase of a DA system, but for new observations to be of value in operational meteorology, they must add positive impact on top of all the other observations already being used, which is much more difficult to demonstrate.Objective verification of precipitation is notoriously difficult, mainly due to the lack of observation data with sufficient density in space and time in most countries.This is improving, mainly due to high-resolution observing techniques such as weather radars providing 2-D precipitation data, but for the majority of the period covered in this review, subjective verification (comparison of precipitation maps to point observations) and pointbased verification have been the main tools.
Table 5 gives a summary of the DA techniques used to assimilate the GNSS tropospheric products.It can be seen that all NWP models used operationally in Europe have the capacity to assimilate ZTD, with the exception of the COSMO model, which is nudging based and uses IWV.
Most of the early studies (i.e.before 2005), showed a positive impact on short-term precipitation forecasts from assimilation of GNSS data.They are often most clearly seen by visual inspection of precipitation patterns, in particular regarding case studies of heavy precipitation events.This is often associated with improvement in the humidity (compared to radiosonde and ground measurements) and of cloud cover.Occasionally negative effects are observed, e.g. a drying of the COSMO model in winter and increased overprediction of light rain (DMI HIRLAM), etc.
It was the UK Met Office and Météo France that first started large-scale operational use of GNSS data in Europe, following the establishment of E-GVAP in 2005.Today both institutions report a positive impact from the use of GNSS data in both regional and global NWP models.Several other institutions now use GNSS delay data operationally, and it is expected that many more will follow in the near future.
Several institutions are now starting to run hourly rapid update cycle (RUC) models with hourly DA.Depending on the cut-off time, this requires faster access to GNSS ZTDs delivered from sub-hourly or RT GNSS processing systems.Such applications can be referred to as "NWP nowcasting".It does not differ greatly from standard NRT NWP, except for the increased demand for higher-resolution, high-frequency observations.In some cases, the cut-off times are reduced so much that many "standard" observations have not arrived at the NMHSs when the RUC DA is run, in which case the RUC model will need to rely more heavily on specific frequent fast access observations, such as from sub-hourly GNSS (de Haan, 2013).

Handling GNSS observation errors, variational
bias control, observation error correlations, and data thinning The amount of ground-based GNSS data is rapidly increasing, which calls for a higher degree of automatisation in estimation of observation errors in NWP vs. GNSS biases.Schemes for VAR bias control (BC) are being introduced and tested in both AROME (Moll et al., 2008)  There is a huge range of NWP resolutions involved; grid box sizes range from ∼ 60 km (ARPEGE on the far side of the Earth viewed from Paris) to 1.5 km in the latest UK Met Office local model.Many NMHSs now run their highresolution models with a ∼ 2.5 km grid box size (AROME, COSMO, HARMONIE), and many will move toward 1 km resolution in the near future.The NWP resolution is important for the optimal spatial density of the observations presented to the NWP DA system.Recently Météo France reduced their thinning distance to 10 km for AROME at a 2.5 km grid box size, meaning effectively no data thinning of interdistances at typical present-day GNSS sites.This increased the number of sites assimilated by a factor of 4, and measurably improved the forecasting skill of the AROME NWP model at Météo France (Moll et al., 2008).More recently, Szintai and Mile (2015) showed that for AROME over Hungary, the effect of humidity observations per observation on the data analysis is larger than for other observations, and largest for GNSS ZTD, which indicates that the DA system will benefit strongly from more of these type of observations and that we are far from saturation as regards the density of GNSS observations in the system.Similar indications have previously been found in case studies, where additional field campaign data were added to the operational GNSS ZTD data.
The quality of GNSS ZTD estimates systematically improve with time, from RT to NRT, and further on to postprocessed data.Within operational meteorology this is a highly uncommon issue, and meteorological databases and DA systems are not built to handle quality updates of the same observation type.In E-GVAP, this is currently handled by naming the solutions differently, e.g.METR vs. METO (sub-hourly vs. hourly).Handling this issue in a more sys-tematic way, transparent to all users and data producers, will likely improve the impact of GNSS tropospheric products in NWP.The importance of this increases with the movement toward rapid-update NWP with short cut-off times, while still running global NWP with a 12 h DA time window.
It is clear from the nature of GNSS processing that GNSS ZTDs will have correlated errors in space and time.However, in NWP DA it is normally assumed that observations do not have correlated errors, which simplifies and speeds up the DA.If one assimilates observations with correlated errors, the data will have a higher influence on DA than they deserve according to their individual observations errors.While some work has been done to estimate the GNSS ZTD error correlations, very little work has been carried out to try to handle the correlations in DA (see TOUGH Deliverables D19 and D22, http://tough.dmi.dk/deliverables/d19.html and http: //tough.dmi.dk/deliverables/d22.html).A simpler method is to reduce error correlations by data thinning, but as mentioned above, indications are that more GNSS ZTDs are beneficial to NWP DA.

Assimilation of SPD and ZTD gradients
Assimilation of ZTD gradients and SPD appears to be a natural next step in the NWP use of GNSS tropospheric observations.
In TOUGH, a 3-D-Var assimilation operator for SPD was implemented in the HIRLAM NWP model (by R. Eresmaa, see TOUGH Deliverable D38, R. Eresmaa, http://tough.dmi.dk/deliverables/d38.html).Several NMHSs worked on enabling slant estimation and a number of impact experiments were carried out.Some of this is documented in Eresmaa et al. (2007) and Järvinen et al. (2007), who tested assimilation of both synthetic and real SPD data.Synthetic data were found to result in reasonable analysis increments.However, the assimilation of real SPDs was found to result in too large analysis increments; however the increments could be reduced by reducing their weighting to correct for error correlations between numerous slants from a single receiver.With the weights reduced, the impact on the analysis of specific humidity corresponds well to that from radiosonde in the same region.Zus et al. (2011) implemented a SPD assimilation capability to the 4-D-Var MM5 model and verified its functionality.Assimilation experiments indicate that the impact on the precipitation forecast is weak, but positive.
Since COST-716 and TOUGH, the typical NWP resolution at operational NMHSs has increased significantly and is likely to benefit from use of asymmetrical GNSS products, such as ZTD gradients and SPD estimates.3.2 GNSS for monitoring severe weather events and operational nowcasting Severe weather events are very diverse, from very cold winters to very hot summers, and the latest suggestion is that they are being observed more and more frequently (Field et al., 2012).One type of severe weather event is intense precipitation, often associated with strong convection potentially leading to flash floods, large economic losses, and risk to life.Over the last few years, several publications discussed the use of GNSS for studying water vapour field evolution, and the production of GNSS-based products tailored for severe weather event monitoring (Table 6).
Several authors have studied the relationship between water vapour fields and the life cycle/intensity of precipitation systems.In general they report (1) a clear relationship between the variations of the GNSS IWV, surface pressure, and precipitation intensity (e.g. in north of Spain during the period 2002-2010 Seco et al., 2012), (2) a preference for frontal systems to develop where the largest amount of water vapour is available, and (3) evidence that water vapour has a predominant role as a precursor for initiation of local convection (van Baelen et al., 2011).These convective systems are associated with typical water vapour configuration characteristics (such as a dry/wet contrast), changing rapidly over time.They can also be initiated due to the orography, causing a large transfer of water vapour, e.g. between plains and mountainous regions (Graham et al., 2012), visible in GNSS IWV, leading to convection and producing thunderstorms.
Various studies have demonstrated the capability of GNSS to monitor small-scale water vapour structures associated with the initiation and development of convective systems: the link between the 2-D moisture field (2D GNSS water vapour maps), the activity of the thunderstorm and lightning is clearly established (de Haan, 2008) and can be used to identify situations of deep convection condition and provide nowcasting warnings (Brenot et al., 2013).The use of GNSS IWV to nowcast lightning activity was also investigated (Mazany et al., 2002), where they demonstrated the capability of GNSS IWV combined with other meteorological parameters to provide a lightning prediction index valid for the next few hours.
In their studies, all these authors obtained clear evidence of the benefits that GNSS can bring to the monitoring of severe weather events.They also often mention that further research and developments are required to improve and standardise their methodology and to be able to allow use in operational nowcasting at NMHSs.The availability of the ultra-fast/real-time GNSS-tailored products (Sect.2.2) will undoubtedly boost this research field in the near future.

GNSS tomography
The potential of dense GNSS networks to derive spatially resolved humidity fields, and the application of tomographic techniques was already emphasised by Ware et al. (2000) and Bevis et al. (1992).Bevis et al. (1992) suggested a setup of small dense GNSS networks specifically designed for meteorological applications.Such networks have so far only been established for specific campaigns.However, nationwide geodetic GNSS networks are continuously expanding and the station density is in some sufficient regions for tomographic applications.Currently, there are active GNSS tomography groups in most European countries which develop and operate tomography systems (Brenot et al., 2014;Benevides et al., 2014;Reverdy et al., 2009;Champollion et al., 2009;Notarpietro et al., 2011;Miidla et al., 2008;Perler et al., 2011;Kačmařík and Rapant, 2012;Rohm, 2013;de Haan, 2008;Saqellari-Likoka and Karathanassi, 2008).
The basic idea of GNSS tomography is to combine a large number of signal delays, which cover most parts of the atmosphere from different directions in order to obtain a spatially resolved field of the atmospheric refractivity or humidity.A single SPD is defined by where the geometric delay S −G = s ds − g ds describes the extra path length due to ray bending.This term is usually neglected in GNSS tomography.The refractivity N is related to the atmospheric quantities pressure, temperature, and humidity, and N (s) is the refractivity along the bent signal path.
To obtain the 3-D refractivity field, a large number of SPDs are required which are combined with the observation vector m.The 3-D refractivity field is then mapped to the state vec-tor x and Eq. ( 5), which relate each single observation to the state x, are combined to a vector operator F .This leads to a non-linear vector equation which needs to be solved for x, i.e. the 3-D refractivity field.
The equation is non-linear as the bent signal paths depend on the unknown refractivity field.This leads to severe problems and a linearised form is derived where the signal paths are assumed to be known and to be independent from x.The liberalisation does not necessarily require the signal paths S to be straight lines but this is usually assumed.To solve the equation numerically, the problem is discretised and the 3-D refractivity filed is defined on the nodes of a spatial grid.
The linearised and discretised form of the integral in Eq. ( 5) leads to a weighted sum of the N j along each signal path, where the weights specify the contribution of each N j to the signal delay i: SPD = m i = j ω j N j .The weights ω j can be combined to a matrix A, which consists of one row i for each observation and one column j for each grid node, and Eq. ( 6) can thus be written as a matrix vector product Only the weights ω j , i.e. the matrix elements a ij , near the signal path contribute to a given SPD, all other matrix elements are filled with zeros.Consequently, the matrix A is a large sparse matrix.In general, A is not a square matrix as the number of observations is not related to the number of grid nodes which define the refractivity field.The inverse matrix A −1 does, therefore, not exist and a least-squares solution does not provide a stable solution because A is a sparse illconditioned matrix.This is typical for ill-posed problems which do in general not have a unique solution and the solutions are usually not stable.Regarding Eq. ( 7), there is in general an infinite number of different vectors x which would lead to the same set of observations m and small variations in the observations m, e.g.observation errors, do lead to completely different solutions x.Tomography generally leads to ill-posed problems and a large variety of solution strategies have been developed.However, most of these strategies are not optimal for GNSS tomography.The GNSS ground networks used today were designed for geodetic applications with rather large interstation distances and a inhomogeneous distribution of stations.The satellite-receiver links which scan the atmosphere vary with the satellite constellation and are highly variable in space and time.The properties of the operator A are, therefore, also time dependent, which puts high demands on the inversion strategy.
The solution of Eq. ( 7) can be separated into four subtasks: Discretisation The spatial discretisation of the atmospheric state defines the number of unknowns.A large number of N j would be optimal to obtain a good approximation of Eq. ( 5) while a minimum number of unknowns would be best for the inversion of Eq. ( 7).The grid dimensions determine the number of SPDs which are covered by the grid and the vertical distribution of nodes defines the vertical resolution, especially of the important boundary layer.

Modelling of the operator A
There is some freedom in the modelling of the operator A, i.e. in how to choose the weights ω j = a ij which relate the N j to the integral N ds.In general, a smooth operator with a small operator error would lead to more stable solutions.

Inversion strategy
The choice of the inversion strategy depends on the properties of the operator A. There is a large variety of standard techniques which can be used to solve Eq. ( 7).
Constraints The problem defined by Eq. ( 7) is globally underdetermined and additional observations or constraints are required to obtain stable solutions.
These four elements are essential for reliable and stable results but a continuously running tomography system is much more extensive.
The atmosphere is a dynamic system and observations made at different times are related to different atmospheric states.In an ideal case the tomography should use only SPD observations from a certain time but the information provided by such a data set is usually insufficient to define the 3-D refractivity field.The spatial information is, therefore, increased by collecting observations over longer periods, i.e. by reducing the temporal resolution.Depending on the sampling rate of the SPD observations, periods between 30 s (Perler et al., 2011) and 3 h were used.Some inversion algorithms require error estimates of the first-guess state and the observations as well.Depending on the data which enter the tomography, efforts must be made to estimate the errors.
Another important factor is the processing of the SPD observations which defines the observation error to a large degree.Many inversion algorithms show a tendency to amplify the observation errors which gives the SPD processing strategy (see Sect. 2.2) a high impact on the performance of the entire tomography system.As it is not possible to separate the impact of the observation error from other influential factors, this will not be further discussed in this section.
The GNSS tomography experiments reported so far are rather different in many respects, and comparing the quality of the results is almost impossible.Furthermore, tomography results in large 3-D fields but validation studies are usually limited, e.g. to a single radiosonde station or some horizontal layers and the relative performance of different approaches remains an open question.The current state of the European GNSS tomography will, therefore, be reviewed regarding the four criteria mentioned above.An overview about process-ing and inversion strategies can also be found in Notarpietro et al. (2011).

Spatial grid
The spatial discretisation, i.e. the choice of a specific grid is an important issue and depends on the application.Large nationwide networks (Troller et al., 2006a;Bender et al., 2011a) require different grid to small limited-area networks.
In most experiments, a regular horizontal grid is used, sometimes extended by a rather coarse outer grid around the region of interest (Troller et al., 2006a;Notarpietro et al., 2011).SPDs at low elevations propagate through large parts of the atmosphere until the top layer of the grid is crossed, e.g.> 50 km at an elevation of 10 • and a model height top of 10 km.Extending grids around small networks is, therefore, a strategy for processing SPDs at low elevations, which provides most of the information regarding the vertical structure.It was found that the horizontal grid spacing should not be much smaller than the mean interstation distance (Brenot et al., 2014;Bender et al., 2011a).The vertical grid spacing is usually chosen to increase with height, thereby covering the important boundary layer with many cells but an equidistant spacing is also possible (Brenot et al., 2014;Bender et al., 2011a).The vertical resolution is limited to several hundred metres, ranging from 200-500 m near the ground, to 1-2 km at higher altitudes.It has also been reported that the absolute location of the grid relative to the GNSS stations has a rather large impact on the results (Notarpietro et al., 2011;Chen and Liu, 2014).

Modelling the operator
The discrete linear operator, i.e. the matrix A, maps the state x on the observations m and the matrix elements a ij specify the contribution of each grid node to a certain SPD.One method of estimating the a ij is provided by Eq. ( 5): the segments j in N ds ≈ N j j can be regarded as the subpaths in each grid cell (= voxel) and N j is the constant refractivity in each voxel.The latter assumption leads to a step function along the signal path and to rather large operator errors, especially for coarse grids.This approximation is used by almost all groups.However, discretisation is not necessarily related to a voxel structure.The quantities at the grid nodes can always be interpolated on any point on the signal path, providing a much smoother profile along the path and leading to much more realistic mapping from the grid to the SPD and the error of the discretisation is reduced (Perler et al., 2011).More recently, Zus et al. (2015) developed operators, which are both fast and accurate, take ray-bending into account, and are thus tailored for variational data assimilation/travel time tomography.

Inversion strategy
Different inversion strategies are available to solve Eq. ( 7) and many of them were successfully applied to GNSS tomography.
The algebraic reconstruction techniques (ART) are a rather simple but fast iterative approach which apply small corrections to an initial state x 0 until a good approximation is obtained (Notarpietro et al., 2011;Bender et al., 2011b).There exists no general criterion for when to stop the iteration, (e.g. when the best approximation is found and oscillating or degrading situations need to be avoided), and in most cases a good initial state is required to obtain stable results.
The pseudo inverse is an attempt to invert Eq. ( 7) which is based on the singular value decomposition (SVD) (Rohm, 2013;Flores et al., 2001).In case of ill-posed problems, the SVD leads to a large number of very small singular values (large condition number of A) which need to be truncated at an arbitrary level.
A weighted, damped, least-squares solution of Eq. ( 7) can be found if an initial state is introduced (Reverdy et al., 2009;Baelen et al., 2011;Brenot et al., 2014;Benevides et al., 2014;Kačmařík and Rapant, 2012).This is equivalent to variational data assimilation where the difference between the observations and the initial state is minimised (de Haan, 2008).A more general approach is the Kalman filter which considers the errors of the initial state and the observations in a consistent way to provide the error of the final state (Perler et al., 2011;Rohm et al., 2014;Champollion et al., 2009;Miidla et al., 2008;Nilsson et al., 2007).The Kalman filter seems to be the best established strategy but it leads to rather large matrices if a large number of observations need to be processed or if the number of grid nodes becomes large.As the inverse of a large matrix needs to be computed, this can lead to rapidly growing processing times.
The inversion strategy might provide an optimal result, but in the case of globally underdetermined problems, the optimum is very often a rather poor refractivity field with lots of artefacts.Up to now, no superior inversion algorithm could be identified and good results were obtained with all of them but, because of insufficient observations, all algorithms produce severe artefacts.

Constraints
Beyond the basic algorithms some constraints are required, as the problem is globally underdetermined, and the solution is non-unique.The coverage of the atmosphere by slant paths is currently rather poor and many grid cells usually do not contain any observations.Therefore, some extra information needs to be provided to obtain stable results and to identify a meteorologically reasonable solution.This can be achieved in many different ways.
Intervoxel constraints are used to distribute the SPD information within the grid and to smooth strong gradients.They are usually realised as filters which replace the quantity at a specific node by a weighted sum over some neighbouring region.Covariance functions with inverse distance weighting (Troller et al., 2006b) or Gaussian filters (Bender et al., 2011b) are then implemented.
Boundary values help to keep the results within a meteorologically reasonable range and to reduce artefacts.The humidity above the tropopause is rather low and the top layer of the tomography grid between 8 and 15 km can be set to zero humidity (Troller et al., 2006a).The surface layer (Bender et al., 2011b) or the all layers below a given height (Troller et al., 2006a) can be constrained by synoptic observations.
The SPD distribution is very inhomogeneous which can lead to unstable results full of artefacts.Pseudo-observations are a way to obtain a more homogeneous SPD distribution within all parts of the grid.Pseudo-observations are a combination of real observations and physical models and can be used to reduce the degree of ill-posedness of the inverse problem.As they do not add any real information about the atmospheric state they usually have a strong smoothing effect.
The deficiencies of GNSS slant observations can be reduced by adding real observations from different sources such as synoptic observations, radiosonde profiles, any kind of profiler data, radio occultation profiles or even GNSS ZWD/IWV data (Champollion et al., 2009).Extra data usually have a stabilising effect on the inversion in the same manner as constraints.
Most groups use combinations of constraints and additional observations to obtain stable and reliable results.Using a carefully tuned tomography system, a continuous series of stable results can be obtained providing important information of the temporal evolution of humidity fields (Brenot et al., 2014;Labbouz et al., 2013;Manning et al., 2012;Baelen et al., 2011).

Prospects of GNSS tomography
First GNSS tomography experiments were reported in the late 1990s (Braun et al., 1999;Foelsche and Kirchengast, 2001), describing the set-up and results of rather small localised GPS networks.Today, after more than 15 years of experience, the networks are enlarged up to nationwide dimensions and numerous inversion strategies have been tested but GNSS tomography is still in an experimental phase.There are two main reasons for the rather slow progress: (1) the number and density of observations is still insufficient compared with the dimensions of the atmosphere and the short correlation length of water vapour and (2) the inversion strategies investigated so far are not adequate to deal with the highly variable problem of GNSS tomography with an almost random distribution of GNSS stations and a fast changing satellite constellation.The degree of ill-posedness is changing with time and the inversion needs to address the variable character of the underlying mathematical problem.
The standard techniques used so far are not optimal in this respect.
The future prospects of GNSS tomography, however, are encouraging.Regarding the development of the number of GNSS sites processed by E-GVAP within the last 10 years and the improvements in processing quality and speed, it can be assumed that the density of GNSS stations and the quality of SPD data become sufficient for GNSS tomography in the near future.Further improvements can be expected when GLONASS, BeiDou, and Galileo tropospheric products become available.Advanced adaptive inversion techniques are currently being developed and need to be applied to GNSS tomography.Real-time humidity fields might be available in the future on a national or even a European scale if efforts are made to densify GNSS networks and to process SPDs in a consistent manner.

GNSS for climate monitoring
Atmospheric water vapour is the most abundant greenhouse gas involved in the climate feedback loop.As the temperature of the Earth's surface and atmosphere increases, so does the moisture-holding capacity of the atmosphere, and atmospheric water vapour is expected to increase in a warmer climate.Both climate models and observations indicate IWV increases by about 7 % per 1 • C increase of temperature (Wentz and Schabel, 2000;Allen and Ingram, 2002;Trenberth et al., 2005;Held and Soden, 2006).Connected to such changes in the water vapour content are changes in the hydrological cycle, i.e. evaporation and precipitation (Bengtsson, 2010).The gradient in evaporation-precipitation increases proportionally to the lower tropospheric water vapour leading to larger differences between dry and wet areas.As a consequence, good knowledge about the regional distribution of the water vapour content of the atmosphere is crucial as it ultimately determines the rate of precipitation.
Traditionally radiosondes and satellite observations have been used for IWV trend analyses.Gaffen et al. (1992); Ross andElliott (1996, 2001) found an upward trend in IWV obtained from radiosonde observations from 1973 to 1995 over North America, except for the north-eastern part of Canada.They also found increases over China and the Pacific Islands.Over the rest of Eurasia, a mixture of positive and negative trends were found, with a tendency for negative trends over eastern Europe and western Russia.Trenberth et al. (2005) found that ocean satellite observations have a positive IWV trend of 0.40 kg m −2 decade −1 from 1988 to 2003.However, despite a large radiosonde network, the temporal resolution is low and differences in calibrations can lead to systematic errors in humidity (Wang and Zhang, 2008;Wang et al., 2013).Some remote sensing methods, observing in the infrared and the optical frequency bands, are limited to clear sky conditions.Other methods, using microwave remote sensing techniques, can be used also during cloudy conditions.On the www.atmos-meas-tech.net/9/5385/2016/Atmos.Meas.Tech., 9, 5385-5406, 2016 other hand, they only provide usable water vapour observations over oceans (Trenberth et al., 2005).GNSS and DORIS can provide IWV observations with high accuracy and complement radiosonde and satellite measurements (Bock et al., 2007(Bock et al., , 2010;;Wang et al., 2007).As the time series of GNSS data grows longer they can also be used for the detection of trends and other systematic effects.
The capability of using GPS data to monitor climate change (e.g. as a linear trends in IWV) has been investigated in only a few studies.Early attempts were made by Gradinarsky et al. (2002), who used data from 1993 to 2002 in Sweden, and found very small linear trends.In their study, the largest trend was seen at Onsala (0.2 kg m −2 year −1 ), on the Swedish west coast, where nearby microwave radiometer data and radiosonde data confirmed the GPS estimates.The study was completed by Nilsson and Elgered (2008); Jarlemark et al. (2010) and Ning and Elgered (2012).These authors found IWV trends in the range from −0.5 to +1.0 kg m −2 decade −1 for the same area, but for different periods.They also studied several sources of uncertainty in the IWV trends due to natural variability of the IWV and errors in the GPS data (antenna phase centre variations, cutoff angle, multipath).Other studies provided IWV trend estimates and comparisons with independent data over specific regions (Morland et al., 2009;Sohn and Cho, 2010) or globally (Jin et al., 2009;Heise et al., 2009;Vey et al., 2010).Most notably, Vey et al. (2010) found good agreement between GPS estimates of seasonal and interannual variations of IWV and NCEP (National Centre for Environmental Prediction) reanalysis, except in the tropics and in Antarctica where the reanalysis underestimated IWV by 40 and 25 %, respectively.Bock et al. (2014), also reported global IWV trends in the range ±2 kg m −2 decade −1 from GNSS and DORIS data, which are in good agreement with ECMWF reanalysis (ERA-Interim) and microwave satellite data (over the oceans).Vey et al. (2009) also investigated the homogeneity of long-term GPS IWV time series and found offsets of up to ±1 kg m −2 due to antenna and radome changes, sudden changes in the number of observations (usually associated with failures in equipment), or changes in the observation cut-off angle used in the GNSS processing.Changes in antenna phase centre or a varying environment affecting the signal multipath will affect the absolute level of the estimated IWV. Ning and Elgered (2012) showed that for some stations it can be advantageous to use elevation cut-off angles as high as 25 • in order to reduce such effects.The idea is that when searching for long-term trends, it can be acceptable that the uncertainty of individual estimates is increased if one at the same time reduces the size of systematic errors.
A recent comparison (Ning et al., 2013) using a regional climate model shows interesting systematic patterns (in both the phase and the amplitude) between coastal and inland sites, where the GNSS results in general confirm the variation seen in the climate model, although differences that call for more detailed studies also exist.
The application of GNSS for climate monitoring is an emerging new field of research.Ongoing reprocessing efforts within EUREF and IGS (see Sect. 2.3) using state-of-the-art models will provide a consistent time series of tropospheric data, taking benefit of more than 20 years of GNSS observations from over 400 stations worldwide.One of the goals of the COST Action ES1206 (GNSS4SWEC) is to assess existing data sets of reprocessed GNSS ZTDs for the period from 1994 to the present, standardise the ZTD to IWV conversion procedure and produce a homogenised data set.This unique data set will (1) enable validation and quantification of systematic biases from a range of instrumentation, (2) improve the knowledge of climatic trends, and (3) benefit both global and regional NWP reanalyses and climate model simulations through assimilation of GPS data or for validation purposes.

Conclusions
The meteorological application of GNSS is now a wellestablished technique in Europe.Currently, the operational NRT tropospheric products, provided through E-GVAP, are used for validation and assimilation in a number of operational NWP models, used in reconstructing 2-D and 3-D water vapour fields and in impact studies of convection.However, there are yet further benefits to be obtained from innovation of the current operational GNSS tropospheric products.At present, most E-GVAP ACs process GPS-only data in a network solution providing hourly updated ZTD-only tropospheric products.Advanced GNSS tropospheric products such as horizontal ZTD gradients, SPD (delays in the direction of each satellite), and 3-D refractivity or humidity fields (using tomographic reconstruction) can now be produced.ZTD gradients and SPD can be derived for receivers independently of their spatial density, while tomography requires input from a dense network in order to function well.
Furthermore, multi-GNSS processing will improve the accuracy of tropospheric products due to the increased number of observations and improved coverage of azimuth and elevation angles.Development and testing of new multi-GNSS products is of particular interest to operational NWP and for the forecasting of severe weather.Short-term, highresolution NWP forecasting or nowcasting models require more detailed humidity observations than are currently available, especially to resolve small-scale phenomena like deep convection.Of particular interest in Europe are regional meteorological extremes such as heavy precipitation events, flash floods and heat waves, which are expected to increase in the future as a result of global warming (Stocker et al., 2013).Relevant areas of research include (1) severe weather forecasting (new GNSS products are required to provide additional information on the spatial heterogeneity and rapid temporal variability of tropospheric humidity), (2) nowcasting (providing rapid updates for the analysis of the atmospheric state requires a transition from NRT GNSS processing as im-plemented in E-GVAP to sub-hourly or real-time (PPP) processing schemes), and (3) multi-GNSS analysis combining data from GPS, GLONASS, Galileo, and BeiDou is expected to provide improved estimates of tropospheric products.
The use of GNSS estimates of the water vapour in climate research has been limited thus far, despite the good consistency between IWV trends derived from GNSS and from other techniques (Ning and Elgered, 2012;Bock et al., 2014).Such results do, however, require homogeneous data sets.Over time, improvements in GNSS processing algorithms have introduced inconsistencies in long-term time series, making climate trend analysis challenging.Ongoing reprocessing efforts within the IGS and EUREF using state-ofthe-art GNSS processing will provide a consistent time series of tropospheric data.These unique data sets will (1) enable validation and quantification of systematic biases from a range of instrumentation, (2) improve the knowledge of climatic trends of atmospheric water vapour, which currently have a large scatter, and (3) benefit both global and regional NWP reanalyses and climate model simulations (e.g.IPCC AR5).
More than 15 years of GNSS meteorology in Europe has already achieved outstanding cooperation, overcoming difficulties such as cross-border data access.Through collaboration between atmospheric and geodetic communities in Europe, it is now feasible to develop next-generation GNSS tropospheric products and applications that can enhance the quality of weather forecasts and climate monitoring.This work is carried out within COST Action ES1206 advanced global navigation satellite systems tropospheric products for monitoring severe weather events and climate (GNSS4SWEC) (2013)(2014)(2015)(2016)(2017).GNSS4SWEC is targeting improved understanding of atmospheric processes that will result in more accurate nowcasting of severe weather.This will lead to improved hazard management, lowering the risk of loss of life and the risk to national infrastructure.GNSS4SWEC will also promote the use of reprocessed longterm GNSS-based tropospheric delay data sets for climate research, with a focus on climate-sensitive regions such as higher latitudes.A correct representation of water vapour in climate models is essential for improvement of global warming and precipitation projections, on which the development of socio-economic response strategies are based.Also, the benefit to NWP and climate modelling in turn will lead to improved satellite-based positioning, navigation, and timing products through improved signal propagation models.In addition, while the production, exploitation, and evaluation of operational GNSS tropospheric products for NWP is a wellestablished research field in northern and western Europe, it is still an emerging R&D field in eastern and south-eastern Europe.GNSS4SWEC is supporting the knowledge transfer and establishment of new ACs, which will benefit the operational E-GVAP service.

Figure 1 .
Figure 1.European GNSS meteorology projects, from 1996 to the present.

Figure 2 .
Figure 2. Number of ground-based GNSS stations delivering NRT ZTDs for operational NWP within E-GVAP.

Figure 3 .
Figure 3. Schematic presentation of individual slant path delays (SPDs) from three GNSS satellites and their mapping to zenith total delay (ZTD).

Table 1 .
Characteristics of the different GNSS tropospheric product type.

Table 3 .
Selected comparisons of GNSS water vapour related studies.The different techniques are further described in the text. 1 zenith wet delay, 2 integrated water vapour, 3 zenith total delay. *

Table 4 .
Comparison of GNSS tropospheric products and NWP models in Europe.Schwitalla et al. (2008)* The different models are further described in the text. 1 integrated water vapour, 2 zenith total delay, 3 slant wet delay.

Table 5 .
Assimilation of GNSS tropospheric products in NWP models.
(Arriola et al., 2016)t al., 2016)NWP models.VAR BC enables use of temporally variable biases in the DA, which is likely a benefit, since errors and biases are often weather dependent.It is also well established that O-B offsets have a seasonal dependence.

Table 6 .
Case studies of severe weather and operational nowcasting products.