Last modified on Feb 27, 2024 15:05

Contributors: Kalev Rannat (Tallinn University of Technology, TUT), Hannes Keernik (TUT), Fabio Madonna (CONSIGLIO NAZIONALE DELLE RICERCHE – ISTITUTO DI METODOLOGIE PER L'ANALISI AMBIENTALE, CNR-IMAA), Emanuele Tramutola (CNR-IMAA), Fabrizio Marra (CNR-IMAA)

Issued by: CNR-IMAA / Fabio Madonna

Date: 03/06/2021

Acronyms

Introduction

This document provides information about GNSS measuring principles, hardware and raw data processing. It is a supplementary document to the Product User Guide (PUG) and Specification for GNSS IPW providing a more extensive overview on the workflow of GNSS IPW product for the CDS.

The fundamental observable of GNSS is the signal propagation time from the satellite to a GNSS receiver. Due to the atmosphere, the signals emitted by the satellite constellations (GPS, GLONASS, Galileo, BeiDou, QZS, etc.) are delayed introducing positioning errors if not handled properly. The tropospheric delay consists of two parts: “dry” delay and “wet” delay. The first one depends on dry gases between a GNSS receiver and the satellite and can be estimated from meteorological parameters analytically. The latter, depending on water vapour amount along signal’s path, is difficult to estimate due to its variability and relatively low contribution (around 10%) to the tropospheric total delay. Nevertheless, the wet component can be calculated using the tropospheric product from GNSS data processing supported by precise near surface meteorological data at each station. Since 1980s, when several studies showed the feasibility of sensing the atmosphere using the by-product of geodetic positioning i.e. the signal delay (e.g. Askne and Nordius, 1987 [1]), GNSS meteorology as a new interdisciplinary research field started to develop [2]. During the last decades, the ability to monitor and interpret GNSS signal propagation delays has proved to be a useful tool in numerous applications. For instance, it is widely used in investigating extreme weather events, detecting changes in climate and improving the initial fields of numerical weather prediction (NWP) models. An overview of the future prospects of the GNSS meteorology can be found in Guerova et al., 2016 [3].

Both geodetic processing and GNSS meteorological data (IPW and ZTD) assimilation methodologies are currently under implementation at several institutes. For example, in Europe, there is a continuous flow of GNSS observations processed in near-real-time from more than 500 stations by more than 10 geodetic analysis centres.

The complexity of the methodology behind GNSS technique requires an extensive knowledge based on different disciplines such as geodesy, meteorology and physics (not to mention notable computational resource required to perform the analysis). Therefore, in this document we provide a coherent description of the GNSS technique’s background (including state-of-the-art in geodesy, satellite navigation, propagation of electromagnetic wave, etc.) although the readers interested in the basics are encouraged to have a look on different sources in the literature review. Comprehensive scientific literature [4, 5, 6] gives an overview about GNSS basics, while GNSS IPW uncertainty budget is discussed in Ning et al., 2016 [7].

The document is organized as follows. The first sections provide an introduction to the GNSS data acquisition with the most common techniques and software (Section 2). Then, a short overview about the signal propagation and tropospheric refraction is given (Section 3). The retrieval of GNSS tropospheric product is described in Section 4. Then, retrieval and usage of ancillary near surface meteorological and site metadata is described together with the retrieval of GNSS IPW product for the CDS, along with its uncertainty (Section 5).

As described in the PUG, only data at their current level of harmonization from the IGS (accessible from CDDIS) are considered for the generation of the IPW product [8]. In addition, IPW derived from the EPN-repro2 dataset, being a reprocessed product and based on data from EUREF Permanent Network (EPN), is made available to the CDS users. Any extension to the GNSS datasets provided through CDS will be considered in future C3S activities. Neither the routinely uploaded IGS tropospheric product nor the historical EPN-repro2 dataset provide any estimation of the IPW and its uncertainty: the algorithm developed in the C3S activities allows the estimation of these quantities and their harmonization applying a consistent methodology (described in Section 5) at all IGS and EPN sites.

Instruments for GNSS and meteorological data acquisition

Principles of GNSS observations

GNSS-estimated IPW is derived from monumented instruments that use GNSS (GPS and other satellite navigation systems) information. The method is based on detecting the signal propagation delays on the path between the transmitter (GNSS satellite) and receiver (GNSS-receiver) caused by atmospheric refraction.
Assuming that the GNSS satellite orbits (positions at each epoch) and the GNSS receiver's antenna locations are known and that the receiver and satellite clocks synchronized (Figure 1), it is possible to calculate the expected propagation time of a signal, knowing the propagation speed in the vacuum of an electromagnetic wave (approximately 3∙10⁸ m/s) and the shortest geometrical path from the GNSS-satellite to the GNSS receiver's antenna. However, the atmospheric refractive index is neither equal to the one in vacuum, nor constant throughout different atmospheric regions. The changes in the refractive index cause two effects. First, the propagation path of the electromagnetic signals gets geometrically bended. Secondly, the velocity of propagation is considerably reduced in the neutral atmosphere. Both effects introduce an excess path delay compared to the time it would take the light to cover the same distance in vacuum. Path delays are estimated for each detected satellite signal at elevation α degrees from horizon. At each measurement epoch, the GNSS receiver gets a number of path delays according to the number of detected GNSS satellites. Then, path delays are converted into Zenith Total Delay (ZTD) by a mapping function (the simplest mapping function is 1/sinα). The final ZTD for a certain measurement epoch and location is a summarized contribution of all path delays detected and mapped into vertical (zenith) direction at the GNSS antenna location. Finally, the ZTD can be decomposed into wet (ZWD) and dry components. The wet component is related to the IPW (description in Section 5 below). A full description of the GNSS data processing from the signal transmission to the calculation of the ZTD including discussion of the accuracy required for precise navigation and geodetic applications can be found in the literature cited in this document.

Figure 1: Illustration of the GNSS-signal propagation in the atmosphere, slant delays (SD) and Zenith Total Delay (ZTD). 

Receivers

In case of reference networks (e.g., IGS and EPN), the usage of receivers and antennas is standardized. Only previously known brands and models as described in the IGS rcvr_ant.tab and IGS20.atx file are accepted with full standing within the IGS network (https://files.igs.org/pub/station/general/ ). There is a wide range of GNSS receiver types that may be used, and the majority consists of a stand-along receiver connected to the internet (either directly or over a gateway). Alternatively, a GNSS receiver may be a PC-card type, e.g., https://www.novatel.com/products/gnss-receivers/oem-receiver-boards/oemv-receivers/oemv-2/

Examples of the GNSS instrumentation and mounting configuration for the antennas at the site are shown in Figure 2.

Figure 2. Geodetic antennas with Dorn-Margolin antenna element and choke-rings (upper left and right) at Onsala and Visby sites (middle left and right), GNSS receiver type SEPT POLARX5TR as used at ONSA and VIS0 sites, the receiver in apparatus rack at site VIS0 (lower left and right) (source: http://www.epncb.oma.be/_networkdata/siteinfo4onestation.php?station=VIS000SWE)

An incomplete list of most widely used receiver manufacturers:
Leica: http://leica-geosystems.com/en-gb/products/gnss-systems
Trimble: http://www.trimble.com/positioning-services/
Javad: https://www.javad.com/
Novatel: https://www.novatel.com/#latestNews
Ashtech: https://www.navtechgps.com/receivers/

GNSS data processing methods

• Network solution (DD). Using Double Differences (DD), the clock errors of both the satellite and receiver are eliminated [5]. A large network is necessary to obtain absolute estimates. Observations of a network of receivers, gathered over a certain time window (e.g. 12 hours) are necessary to determine the position of a receiver accurately. The determination is performed using GNSS processing software, which estimates the position of the receivers in the network and, simultaneously, the atmospheric correction or atmospheric delay.

• Precise Point Positioning (PPP). For this method, the orbits and satellite clocks are estimated using a separate scheme and then used as a priori information to estimate the position of the receiver and atmospheric term [9]. This method requires very accurate and stable satellite information but has the advantage of being completely scalable with respect to the number of GNSS sites in the processing scheme. Both methods give similar quality results if everything is done in a correct and consistent way. The results (in GNSS IPW context the ZTD and its 1σ errors) cannot be classified as “worse” or “better” based on information about the data processing method. For most applications, it is indifferent by which method or software the GNSS data is processed. However, it may be useful for the data analyst to know which method and method-specific constraints were used for more specialized analysis.

Both methods give similar quality results if everything is done in a correct and consistent way. The results (in GNSS IPW context the ZTD and its 1σ errors) cannot be classified as "worse" or "better" based on information about the data processing method. For most applications, it is indifferent by which method or software the GNSS data is processed. However, it may be useful for the data analyst to know which method and method-specific constraints were used for more specialized analysis.

Software for GNSS data processing

The three most widely used geodetic software in scientific communities are as follows:

• BERNESE (http://www.bernese.unibe.ch/ )
• GAMIT/GLOBK (http://www-gpsg.mit.edu/~simon/gtgk/ )
• GIPSY/OASIS (https://gipsy-oasis.jpl.nasa.gov/ )

However, there are far more applications doing the same or similar processing. For example:

• RTKlib - An Open-Source Program Package for GNSS Positioning (https://rtklib.com/ , by Univ. Tokio)

On-line post-processing facilities like:

• AUSPOS (http://www.ga.gov.au/bin/gps.pl )
• Canadian Geodetic Survey CSRS-PPP on-line service (https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php ).

Or, in-house developed solutions, non-commercial, but not open, for example:

• EPOS used by Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ.

Instruments for surface meteorological data acquisition

Ideally, reference quality meteorological sensors should be installed at the GNSS site, as close as possible to the same position and height of the GNSS antenna. In practice, meteorological parameters used for converting ZWD to IWV can come from a variety of sources depending upon availability (in order of preference):

• Co-located reference quality meteorological instruments
• Co-located baseline meteorological instruments e.g. co-located AWS (Vaisala, Paroscientific etc.)
• Neighbouring meteorological site data, adjusting near-surface pressure to the height of the GNSS antenna
• Triangulated/interpolated data from three nearby meteorological surface sites, adjusting near-surface pressure to the height of the GNSS antenna
• NWP models or reanalysis data.

The uncertainty associated with the surface meteorological data must be quantified and accounted for.
Longer distances (30-50 km) between GNSS sensor and meteorological instruments make it difficult to reliably approximate surface meteorological data to the GNSS antenna's geodetic position introducing an additional uncertainty.

GNSS signal propagation and retrieval of ZHD

GNSS signals propagating from GPS satellites to the ground-based receivers are significantly affected by the atmospheric refraction. The upper part of the atmosphere, the ionosphere, consists of charged particles and has a dispersive effect on GNSS signal propagation. The magnitude of the delay depends on the frequency of the signal, typically L1=1575.42 MHz, L2=1227.60 MHz for the GPS. Geodetic grade GNSS receivers process both L1 and L2 signals. The measurements taken at two frequencies give a possibility to use ionosphere-free combination of the signals [4, 5, 6] that allows to neglect the ionospheric effects – the dominant part of the atmospheric refraction. In the following we assume the ionospheric refraction eliminated and pay attention to the tropospheric effects only.

Unlike the ionosphere, the troposphere is a neutral gas (non-dispersive medium) and its refraction cannot be eliminated by multi-frequency techniques. A detailed description of radio signal propagation in the atmosphere, atmospheric refraction and refraction indexes can be found in literature [4, 5, 6, 10 chapter 2.3 “Signal propagation”].

Tropospheric delay can be defined as a delay of the signal propagating through the troposphere compared to a signal propagating with a speed of light in the vacuum. The medium is non-dispersive and the signals are equally delayed. This delay is a function of the tropospheric refractive index [10]. The excess path length (compared to the Euclidean distance between the satellite and GNSS receiver’s antenna) is caused by two effects – the propagation delay and the delay due to the ray bending. The value of the delay depends mainly on the elevation angle the ray arrives to the atmosphere. The ray bending, which is the difference between the real signal path and the straight line between a transmitter and receiver, can be neglected due to the fact that it only accounts for about 3 cm delay at an elevation 10°. Therefore, the tropospheric delay can be written as:

with integration along the geometric path of the signal and the refractivity defined as:

where n denotes the refractive index:

with c denoting the speed of light in the vacuum and v indicates the speed of the GNSS signal along the real ray path.

The refractivity N is mainly determined by three meteorological parameters (atmospheric pressure, temperature and the partial pressure of water vapour) and has been estimated empirically (e.g. [11]). It can be found from different publications as:

where we can distinguish the dry (first term) and wet components (second and third term) of the refractivity. In this equation k₁, k₂ and k₃ denote empirical constants with values k₁ = 0.776 K∙Pa^-1, k₂ = 0.72 K∙Pa^-1 and k₃ = 0.00375∙10⁵ K²∙Pa^-1; p is the atmo­spheric pressure in Pa, T is temperature in K and e is the partial pressure of water vapor.

The constants in Table 1 have been measured experimentally and published by different authors. The numeric values depend on atmospheric composition [12, 13]. Knowing that the atmosphere changes, it must be accepted that the values of these constants change in time (however, the changes are not large, but should be considered for processing historical data). The values of coefficients used for the generation of the CDS datasets can be found in Section 5.

Table 1. Physical constants k₁, k₂, and k₃ found and justified by different authors. Table modified from [14]

As the hydrostatic equation is only valid for the pressure, not for the partial pressure of dry gases, a more useful expression for N is given by Davis et al., [15]:

where ρ denotes the total density of dry gases and water vapour, R_d is the specific gas constant of the dry constituent (R_d = R/M_d, R = 8314.4621 ± 0.0075 J/(kmol·K)), M_d = 28.9644 ± 0.0014 kg/kmol and Z_w^-1 is the inverse compressibility of wet air accounting for the correction of non-ideal gas behavior. The coefficient k’₂ can be determined by using the equation:

where M_w = 18.0152 kg/mol denotes the molar mass of water vapour [15].

The first term in Eq. (5) depends only on the total density and not on the wet/dry mixing ratio. It can be integrated by applying the condition that hydrostatic equilibrium is satisfied [15]:

where g denotes the acceleration due to gravity and ρ is the total mass density at height h, and p is the total pressure. The integration of the first term in N gives the dry (hydrostatic) component of the total tropospheric delay:

where P₀ is the total ground pressure in hPa and g_m denotes the gravity acceleration at the mass centre of a vertical column of the atmosphere by Saastamoinen [16]:

where

where λ and H denote the site latitude in degrees and the height above the geoid in m, respectively.

The real measurements (the signal delays estimated from GNSS-data) are slant delays (like used in tomographic applications). For IPW solution the “slanted delays” are mapped to zenith by a mapping function assuming that the atmosphere is isotropic.

The slant tropospheric propagation delay can be expressed as [15, 17]:

where θ is the elevation angle of the satellite, \( \Delta_{d}^{trop} \)  denotes the zenith dry (hydrostatic) delay (ZHD), and \( \Delta_{w}^{trop} \)  the wet zenith wet delay (ZWD). DryMap is the "mapping function" for the hydrostatic delay (see below) and WetMap is the mapping function for the wet delay.

A mapping function is a mathematical model for the elevation dependence of the respective delays. The mapping functions (for both the dry and the wet terms) are approximately equal to the cosecant of elevation. This approximation assumes a horizontally layered atmosphere and ignores the earth’s curvature. Marini [18] proposed a continuous fraction that describes the geometrical approximation of the elevation dependence of the delay. Different authors have improved the basic Marini mapping functions. The main practical difference between the mapping functions is their elevation angle range for which they are valid with expected accuracy.

Usually, the GNSS data processing software allows choosing between different mapping functions. For example, for meteorological studies, the Global Mapping Function (GMF) developed by Boehm et al., [19] is based on fitting numerical weather model (NWM) data over 20 years. A more accurate reconstruction of the NWM data can be obtained by interpolating hydrostatic and wet mapping function coefficients as a function of time and location from the global grid files compiled by the Vienna group [20], known as a Vienna Mapping Function. There exist also widely used Niell Mapping Functions [21, 22]. The choice between mapping functions is based on user considerations.

The equation for ZHD (the dry component of the tropospheric delay in units of length) can be written as [23]:

This is the practical equation for calculating ZHD at the site, where the atmospheric pressure, site height above the geoid and site latitude are known from the site’s metadata.

The total tropospheric delay (in units of length) can be obtained similarly to the equation for refractivity (Eq. 5) which consists of dry (hydrostatic) and wet components:

Retrieval of ZTD

The first element in GNSS IPW observations’ processing chain (Figure 3) is a geodetic grade GNSS receiver providing L1 and L2 signal phase pseudoranges. A detailed description of specification of GNSS signals and raw data processing technology starting from geodesy, the basics of navigation and signal processing can be found in the literature [4, 5, 6, 24]. A lot of additional theoretical and software-specific information is available in software manuals [17, 25, 26].

Figure 3. Typical GNSS IPW measurements’ main processing chain, where p_s denotes surface pressure, T_s the surface temperature, and Tm the mean temperature of the atmosphere.

L1, L2 pseudoranges (the second item in Fig. 3) are used as input to a Forward Model (geodetic data processing software). The Forward Model (the third item in Fig. 3) can be any of the known geodetic software (see Section 2.4).

In general, the main scientific GNSS packages used for GNSS data processing agree at 5–6 mm ZTD level (~ 1 kg/m² IWV) for instantaneous estimates, or about ~ 2 mm ZTD (0.3 kg/m² IWV) for weekly estimates, and about 0.15–0.30 kg/m² per decade for IPW trends [27].

For describing the origin of GNSS final troposphere product, we start from the observation equations [28]. The ionospheric-free combinations of dual-frequency GNSS pseudorange (P) and carrier-phase observations (Φ) are related to the user position, clock, troposphere and ambiguity parameters according to the following simplified observation equations:

l_P  is the ionosphere-free combination of P1 and P2 pseudoranges (2.546P₁-1.546P₂),

l_Φ is the ionosphere-free combination of L1 and L2 carrier-phases (2.546λ₁Φ₁-1.546λ₂Φ₂),

dT is the station receiver clock offset from the GPS time,

dt is the satellite clock offset from the GPS time,

c is the vacuum speed of light,

T_r is the signal path delay due to the neutral-atmosphere (primarily the troposphere),

N is the non-integer ambiguity of the carrier-phase ionosphere-free combination,

λ₁, λ₂, λ are the of the carrier- phase L1, L2 and L3-combined (10.7cm) wavelengths, respectively,

ε_P , ε_Φ are the relevant measurement noise components, including multipath.

ρ denotes the geometrical range computed by iteration from the satellite position (Xs, Ys, Zs) at the transmission epoch t and the station position (x, y, z) at the reception epoch T = t +ρ/c,

Using the IGS orbit/clock products, the satellite clocks (dt) in Eqs. (14) and (15) can be fixed (considered known) and thus can be removed. Expressing the tropospheric path delay (T_r) as a product of the zenith path delay (zpd) and mapping function (M), which relates slant path delay to zpd, it is possible to obtain the point positioning mathematical model functions for pseudorange and phase observations:

Linearization of the observation Eqs. (17) and (18) around the a-priori parameter values and observations (X⁰, ℓ) in the matrix form becomes:

where A is the design matrix, δ is the vector of corrections to the unknown parameters X, W = f(X⁰, ℓ) is the misclosure vector and V is the vector of residuals. The partial derivatives of the observation equations with respect to X, which in the case of Precise Point Positioning (PPP) consist of four types of parameters: station position (x, y, z), receiver clock (dT), wet troposphere zenith path delay (zpd_w) and (non-integer) GNSS signal carrier-phase ambiguities (N), forming the design matrix A:

with: 

The least squares solution with a-priori weighted parameter constraints (P_X⁰) is given by:

where P_l is the observation weight matrix. Solution of Eq. (20) consists of position coordinates (x, y, z), receiver clock (dT), carrier phase ambiguities (N) and tropospheric path delay (zpd_w). As the Eqs. 14 to 21, given by Kouba 2009, were written for PPP, the final product in our notation is ZWD. In case of IGS, the final troposphere product is ZTD instead (atmospheric delays, the fourth item on Fig. 3). ZTD and ZWD are related with each other using the Eq. 13, where ZHD is estimated from meteorological data (e.g., by Saastamoinen model).

Retrieval of IPW

The GNSS IPW datasets provided through the CDS are based on using tropospheric products delivered by IGS data ACs and EPN-repro2. This section describes the data processing schema for retrieval of GNSS IPW for both IGS and EPN tropospheric product. A simplified description of the data flow and data processing steps for retrieval of GNSS IPW from GNSS tropospheric products is depicted in Fig. 4.

Retrieval of GNSS IPW can be described by 3 main steps:

• Pre-processing of GNSS troposphere product (steps 1–4).
• Pre-processing of ancillary data (steps 5–7).
• Calculating IPW and its uncertainty (steps 8–10).

As a first step, the data – ZTD with its uncertainty – is downloaded from IGS and EPN data repositories (more information about data availability and data access can be found from GNSS PUG). The site metadata for both networks are obtained from SEMISYS [29]. Due to the fact that the site latitude and longitude coordinates are given with height from WGS84 ellipsoid, the height above the mean sea level (AMSL) used in IPW retrieval must be calculated first (step 4, Fig. 4).

During the step 5 (Fig. 4) supporting meteorological data from ERA5 is downloaded and the near-ground pressure, p₀, is interpolated to the site latitude, λ. After this, ZHD is calculated (step 6 on Fig. 4, Eq. 12). According to Eq. (13) ZWD is simply calculated by subtracting ZHD from ZTD.

Relating ZWD to IPW (step 8, Fig. 4) can be done using the equation detailed in [24] (Bevis et al., 1994):

For estimating the IPW we first need to find values for the ZWD (Eq. 22) and the conversion factor Π (step 7, Fig. 4) calculated by:

where

where ρ_w is the density of liquid water, R_w is the specific gas constant for water vapour, T_m is the water-vapour-weighted mean temperature of the atmosphere, M_w and M_D are the molar masses of water vapour and dry air, respectively, and k₁, k₂, and k₃ are physical constants from the widely-used formula for refractivity [30, 31]. The value for T_m is calculated from ERA5 data, namely using vertical profile of water vapour pressure (p_w) and the air temperature (T) along the integration path (S):

T_m is interpolated to the location of GNSS site with bilinear interpolation from ERA5 (step 5 Fig. 4) by Eqs. 26–28. Given four points q₁₁=(x₁,y₁), q₁₂=(x₁,y₂), q₂₁=(x₂,y₁), q₂₂=(x₂,y₂),the first linear interpolation is applied in the x-direction:

Then, interpolation in the y-direction is carried out to obtain the desired estimate.

In case of estimating the surface pressure at the GNSS’ altitude, p₂, pressure conversion through the following hypsometric equation at all four grid points before applying bilinear interpolation is needed:

where p₁ is the ERA5 surface pressure at the grid point investigated in hPa, the difference Z₂–Z₁ is referred to as the thickness between pressure levels p₁ and p₂ in meters and H is the scale height defined by:

where g₀ is earth-surface gravitational acceleration (9.80665 m s^-2), R_d is specific gas constant for dry air (287.05 m² s² K^-1) and virtual temperature, T_v, can be found using:

Here e is the water vapour pressure in hPa, which is given according to the Magnus type equation [32]:

where T denotes temperature (K) and RH represents relative humidity. It must be clarified that ERA5 pressure levels are grid-box indicative estimates and – although not equivalent to the point altitude due to smoothing of the topography – are assumed to be sufficiently accurate. This assumption allows us to calculate pressure at GNSS stations using a simple conversion and a bilinear horizontal interpolation which has sufficient accuracy for GNSS applications (see the papers cited and results presented in PUG).

Since the time resolution of both, ERA5 data and IGS ZTDs, is one hour there is no need to interpolate T_m and p_s in time. On the other hand, this is done in case of EPN-repro2 as its data are provided at half-hours.

For calculating the GNSS IPW uncertainty (step 9, Fig. 4), the approach chosen for the CDS is based on the GRUAN GNSS data processing [7], and accounts for all of the uncertainty sources due to the measurement systems. This approach estimates the total uncertainty of GNSS IPW using the error propagation formula including all the uncertainties of the input variables of the IPW retrieval [33]:

where f(V₁, …, V_M) is the functional relationship between the GNSS IPW and input variables, and σ_i is the related 1-sigma uncertainty. written as:

where the 2^nd and 3^rd component below the square-root are to the uncertainty of the ZHD (Eq. 34).

Uncertainty of the conversion factor Π (Eq. 23) can be written as follows:

where σ_PW, σ_ZTD, σ_p0, σ_c, σ_Π, σ_Tm, σ_k3, and σ_k2’ are the one-sigma uncertainties of IPW, ZTD, p₀, the constant in retrieving of ZHD, the conversion coefficient Π, T_m, and the constants from widely used formula for refractivity k₂^’ = 22.1 K/hPa and k₃ = 373900 K²/hPa respectively, and ρ_w = 1000.0 kg/m³ and R_w = 461.522 J/K∙g are the density of water and the specific gas constant for water vapour. The other constants are σ_p0 = 0.6 hPa (for ERA5), σ_c = 0.0015 (dimensionless), σ_Tm = 1.5 K (statistically determined for ERA5 that is the source of meteorological input data), σ_k3 = 1200 K²/hPa, and σ_k2’ = 2.2 K/hPa.

After the values for IPW and its uncertainty are calculated, the values obtained are ingested to the database in the CDS (step 10, Fig. 4).

Summary

Estimating atmospheric water vapour content based on radio signal propagation through the atmosphere has many advantages. Most notably, it allows monitoring of atmospheric humidity with very high spatial and temporal resolution with reasonable costs and provides measurements that are not affected by rain and clouds.
On the other hand, the technique itself has several requirements and its data process includes many steps. In short, after giving a short overview about the theoretical background of GNSS data processing the four main requirements for GNSS IPW retrieval can be given:

• availability of GNSS tropospheric product,
• availability of GNSS site metadata,
• availability of ancillary meteorological data,
• information about the uncertainties.

The CDS users benefit from the current C3S activity by accessing data that are a combined product derived from several sources: (1) IGS and EPN provide tropospheric GNSS data; (2) SEMISYS provides metadata for all stations within these networks and (3) ancillary meteorological data at GNSS stations is estimated by using ERA5. In addition to applying consistent methodology for estimating IPW based on this information as well as calculating its uncertainty according to the latest approach, the input data used in this process are made available. In such way, the CDS database offers all the data needed for complete traceability of IPW and its uncertainty for more than 700 stations worldwide.

References

[1]    Askne J, Nordius H. Estimation of tropospheric delay for microwaves from surface weather data. Radio Science. 1987; 22:379-386.
[2]    Bevis, M., Businger, S., Herring, T. A., Rocken, C., Anthes, R. A., Ware, R. H., GPS meteorology: Remote sensing of atmospheric water vapor using the Global Positioning System, J. Geophys. Res., 97, 15787–15801, 1992.
[3]    Guerova, G., Jones, J., Dousa, J., Dick, G., Haan, S., Pottiaux, E., Bock, O., Pacione, R., Elgered, G., Vedel, H., Bender, M. Review of the state-of-the-art and future prospects of the ground-based GNSS meteorology in Europe, Atmos. Meas. Tech. Discuss. 1-34, 2016, 10.5194/amt-2016-125.
[4]    Borre, K., Strang, G., Linear Algebra, Geodesy and GPS, Wellesley-Cambridge Press, 1997, ISBN 10: 0961408863, 638 pp.
[5]    Hofmann-Wellenhof, B., Lichtenegger, H., and Collins, J., Global Positioning System, Theory and Practice, 5th Ed., Springer, Wien, 2001, ISBN 3-211-83534-2, 382 pp. 
[6]    Teunissen, P.J.G., Montenbruck, O. (Eds), Springer Handbook of Global Navigation Satellite Systems, Springer International Publishing, 2017, DOI: 10.1007/978-3-319-42928-1.
[7]    Ning, T. et al., The uncertainty of the atmospheric integrated water vapour estimated from GNSS observations, Atmos. Meas. Tech., 9, 79–92, 2016, www.atmos-meas-tech.net/9/79/2016/, doi:10.5194/amt-9-79-2016. 
[8]    Noll, C., The Crustal Dynamics Data Information System: A resource to support scientific analysis using space geodesy, Advances in Space Research, 45(12), 1421-1440, 2010, ISSN 0273-1177, DOI: 10.1016/j.asr.2010.01.018. 
[9]    Zumberge, J. F., Heflin, M. B., Jefferson, D. C., Watkins, M. M., Webb, F. H., Precise point positioning for the efficient and robust analysis of GPS data from large networks, J. Geophys. Res., 102 (B3), 5005-5017, doi:10.1029/96JB03860, 1997. 
[10]    Sansó, F., Satellite Geodesy, 2nd Edition G. Seeber, de Gruyter, Geophysical Journal International, 157:1, 589 pp, 2003, ISBN 3-11-017549-5, https://doi.org/10.1111/j.1365-246X.2004.02236.x. 
[11]    Essen, L., Froome, K. D., The Refractive Indices and Dielectric Constants of Air and its Principal Constituents at 24,000 Mc/s, Proceedings of the Physical Society, Section B, 64:10, pp. 862-875, 1951, DOI: 10.1088/0370-1301/64/10/303.
[12]    Jean M. Rüeger, Refractive Index Formulae for Radio Waves, In Proceedings of FIG XXII International Congress Washington, D.C. USA, April 19-26, 2002, pp. 1–13.
[13]    P. Poli, et al., Meteorological Applications of GPS Signals with Ground and Space-Based Receivers, Inside GNSS, 2008, pp. 30–39. 
[14]    Bevis, M., S. Businger, S. Chiswell, T. A. Herring, R. A. Anthes, C. Rocken, and R. H. Ware, GNSS Meteorology: Mapping Zenith Wet Delays onto Precipitable Water. J. Appl. Meteor., 33, 379–386, 1994, https://doi.org/10.1175/1520-0450(1994)033<0379:GMMZWD>2.0.CO;2. 
[15]    Davis, J. L., Herring, T. A., Shapiro, I. I., Rogers, A. E. E., and Elgered, G.: Geodesy by radio interferometry: Effects of atmospheric modelling errors on estimates of baseline length, RadioSci., 20, 1593–1607, doi:10.1029/RS020i006p01593, 1985. 
[16]    Saastamoinen, J., Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites, in The Use of Artifical Satellites for Geodesy, Geophys. Monogr. Ser, vol. 15, edited by S. W. Henriksen et al., pp. 247–251, AGU, Washington, D. C., 1972. 
[17]    Herring, T. A., King, R. W., McClusky, S. C., GAMIT Reference Guide, Rel. 10.3, Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, 2006. 
[18]    Marini, J. W., Correction of satellite tracking data for an arbitrary tropospheric profile, Radio Sci 7:223–231, 1972, doi:10.1029/RS007i002p00223. 
[19]    Boehm J., Niell, A., Tregoning, P., Schuh, H., The Global Mapping Function (GMF): A new empirical mapping function based on data from numerical weather model data. J. Geophys. Res., 33, L07304, 2006b, DOI:10.129/2005GL025546. 
[20]    Boehm, J., Werl, B., and Schuh, H., Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range Weather Forecasts operational analysis data, J. Geophys. Res., 111, B02406, 2006a, doi:10.1029/2005JB003629.
[21]    Niell, A., Global mapping functions for the atmosphere delay at radio wavelengths, J. Geophys. Res., 101, 3227–3246, 1996.
[22]    Niell, A., Improved atmospheric mapping functions for VLBI and GPS, Earth Planets Space, 52, 699–702, doi:10.1029/95JB03048, 2000. 
[23]    Elgered, G., Davis, J.L., Herring, T.A., & Shapiro, I.I., Geodesy by radio interferometry: Water vapor radiometry for estimation of the wet delay,  J. Geophys. Res., 96, 6541-6555, 1991, https://doi.org/10.1029/90JB00834. 
[24]    Parkinson, B.W., Spilker, J.J., Global Positioning System, Theory and Applications, Vol.1, Progress in Astronautics and Aeronautics, American Institute of Aeronautics and Astronautics, Washington, DC, 1996. 
[25]    Dach, R., Hugentobler, U., Fridez, P., Meindl, M. (Editors), User manual of the Bernese GPS Software Version 5.0, AIUB, 2007.
[26]    Dach, R., Lutz, S., Walser, P., Fridez, P. (Editors), User manual of the Bernese GPS SoftwareVersion 5.2, AIUB, 2015. 
[27]    Jones, J. et al. (eds.), Advanced GNSS Tropospheric Products for Monitoring Severe Weather Events and Climate, Springer 2020, https://doi.org/10.1007/978-3-030-13901-8. 
[28]    Kouba, J., A Guide to Using International GNSS Service (IGS) Products, 2015, https://files.igs.org/pub/resource/pubs/UsingIGSProductsVer21_cor.pdf (last checked 27 April 2023). 
[29]    Bradke, Markus, SEMISYS - Sensor Meta Information System. GFZ Data Services, 2020, https://doi.org/10.5880/GFZ.1.1.2020.005. 
[30]    Smith, E. K. and S. Weintraub, The constants in the equation for atmospheric refractive index at radio frequencies, Proceedings of the IRE, 41(8), 1035–1037, doi:10.1109/JRPROC.1953. 274297, 1953, https://doi.org/10.1109/JRPROC.1953.274297.
[31]    Boudouris, G., On the index of refraction of air, the absorption and dispersion of centimeter waves by gasses, Journal of research of the National Bureau of Standards, Sect. A, Physics and chemistry, 67D, 631–684, 1963. 
[32]    Alduchov, O. A., and R. E. Eskridge, Improved Magnus Form Approximation of Saturation Vapor Pressure, J. Appl. Meteor., 35, 601–609, 1996, https://doi.org/10.1175/1520-0450(1996)035<0601:IMFAOS>2.0.CO;2. 
[33]    Immler, F., J. Dykema, T. Gardiner, D. Whiteman, P. Thorne, and H. Vömel, Reference quality upper-air measurements: Guidance for developing GRUAN data products, Atmos. Meas. Tech., 3(5), 1217–1231, 2010, https://doi.org/10.5194/amt-3-1217-2010. 

Related articles

• Page:
  Nordic Gridded Climate Dataset (NGCD): Product User Guide (PUG)
• Page:
  E-OBS daily gridded observations for Europe from 1950 to present: Product user guide
• Page:
  Climate indices derived from NGCD: Product User Guide
• Page:
  Nordic Gridded Climate Dataset (NGCD): Known issues
• Page:
  Nordic Gridded Climate Dataset (NGCD): Algorithm Theoretical Basis Document (ATBD)