Beyond Mapping Functions and Gradients

Abstract

Mapping functions and gradients in GNSS and VLBI applications were introduced in the sixties and seventies to model the microwave propagation delays in the troposphere, and they were proven to be the perfect tools for these applications. In this work, we revisit the physical and mathematical basis of these tools in the context of meteorology and climate applications and propose an alternative approach for the wet delay part. This alternative approach is based on perturbation theory, where the base case is an exponential decay of the wet refractivity with altitude. The perturbation is modeled as a set of orthogonal functions in space and time, with the ability to separate eddy-scale variations of the wet refractivity.

1. Introduction

The effect of the Earth atmosphere on the propagation of light was noticed just after the invention of the telescope by Galileo Galilei, and tables of atmospheric refraction (bending of ray lights) were already available in the XVII century. After the advent of VLBI observations in the fifties and the launch of the first Earth satellite in the sixties, the modeling of the time delays caused by the neutral atmosphere became a necessity.

The current mathematical structure of the modeling of the propagation time delays, used in almost all GNSS software is given by [1].

$$\begin{array}{cc}\mathit{\delta L}\left(e_0\right)& =m_h\left(e_0\right)\left[L_z^h+\mathit{\text{cotg}}\left(e_0\right)(G_N^{h}cos\left(\varnothing \right)+G_E^{h}sin\left(\varnothing \right)\right]\\ & +m_w\left(e_0\right)\left[L_z^w+\mathit{\text{cotg}}\left(e_0\right)(G_N^{w}cos\left(\varnothing \right)+G_E^{w}sin\left(\varnothing \right)\right]\end{array}$$E1

where $\mathit{\delta L}$ is the slant (extra) delay with respect to propagation in vacuum along the bended ray, $L_z^h$ and, $L_z^w$ are the hydrostatic and wet zenith delays, $e_0$ and $\varnothing$ are the satellite elevation and azimuth angles as seen from the station, respectively; $G_N^h$ and $G_N^w$, $G_E^h$ and $G_E^w$ are the north and east components of the hydrostatic and wet delays gradients; $m_h$ and $m_w$ are the hydrostatic and wet mapping functions. Eq. (1) is used in precise GNSS processing software through the modeling of the phase signal [2, 3, 4].

The mapping functions “map” the so-called slant neutral atmosphere (extra) delay $\mathit{\delta L}$ (i.e. the delay along the bended ray from the observer to the emitter) to two “zenithal delays”, named hydrostatic delay (very often improperly called “dry” delay) $L_z^h$ and wet delay $L_z^w$, essentially caused by the water vapor. The mapping functions are usually written [5] in the form of continuous fractions, that were introduced by Marini [6, 7] and normalized by Herring [8] in the form

$$m\left(e_0\right)=\frac{1+\frac{a}{1+\frac{b}{1+c}}}{sin\left(e_0\right)+\frac{a}{sin\left(e_0\right)+\frac{b}{sin\left(e_0\right)+c}}}$$E2

Other (simplified) forms of the mapping functions can be found in the literature [9], but the mainstream form is always Eq. (2). The gradients themselves, noted by the upper case letter $G$ in Eq. (1) were introduced to compensate azimuthal anisotropic effects [10, 11, 12].

For the last thirty years, the improvements on these formulas mainly focused on better and better determinations of the coefficients a, b, c, by comparisons of these formulas with ray tracing. The literature acknowledges as “best” models the VMF1 and VMF3 families, with some seasonally adjusted coefficients constrained from ray-tracing results with respect to Numerical Weather Models (NWM) [13, 14].

The role of the water vapor in the neutral delay is important, as it can be up to 20% (about 45 cm of the zenithal delay $L_z^w$), with respect to the total zenithal delay $L_z^h+L_z^w$ (about 2.3 m). The other gases, including carbon dioxide, have a negligible role in the neutral delay [9, 15], thus cannot be detected through GNSS processing.

Water is present in its three phases on Earth atmosphere, hydrosphere and continents: solid, liquid and water vapor, with important latent heats between phases. Water vapor in the atmosphere has large sources (evaporation, evapotranspiration) and sinks (rain, snow). Water vapor is also the most important greenhouse gas (beyond carbon dioxide) and the driver of cloud coverage. To describe the water cycle [16] is therefore of the uttermost importance, as evidenced by the so-called Energy Balance models [17, 18] that can be written

$$\begin{array}{c}{C}_1{S}_0\left(1-\alpha \right)=C_2\frac{\mathit{dT}}{\mathit{dt}}+C_3{T}^4\left(1-\beta \right)\end{array}$$E3

Where $S_0$ is the solar constant (1360 W/m²), T is the mean temperature on the Earth surface in Kelvin, t is the time. $C_1$, $C_2$ and $C_3$ are constants.

The coefficients $\alpha$ and $\beta$ are albedos, respectively in the visible and infrared wavelengths, both mainly driven by the water vapor contents of the Earth atmosphere [19]. The coefficients $\alpha$ reflects the cloud coverage, typically today at the 30% level, and the coefficient $\beta$ is an infrared albedo, keeping our planet warm at around 15 °C. Without the greenhouse gases, our planet will be at a freezing mean temperature of −18 °C. They have antagonist effects, an increase of $\alpha$ means a cooling of Earth surface, and an increase of $\beta$means a warming, with a lot of intricacies between the positive and negative feedbacks related to the water vapor cycle of the climate models [20]. The ultimate goal of global long-term climate models [21] is to predict which effect will prevail (this is the dT/dt term in the right side of Eq. (3)).

The study [22] highlights the difficulty of measuring atmospheric water vapor with sufficient spatial and temporal resolution, and with sufficient accuracy, to provide observational constraints. GNSS processing is not the only source of water vapor data in the atmosphere. Remote sensing by satellites is the main provider [23], but the resolution of their data sets is limited by the distance between the satellites and the Earth and their orbital cycles. Besides, satellites are expensive. GNSS receivers, even precise ones, are a lot cheaper, and can provide long-term time series with high temporal resolution. Other ground-based instruments are mainly lidars [24], photometers [25], and water vapor radiometers [26]. The only source providing in situ meteorological data are the radiosondes [27], launched twice per day in a limited amount of worldwide sites. Many studies have been devoted to the causal relationship between water vapor and rain [28, 29], including extreme events [30].

It is therefore important to separate the water vapor modeling coefficients $L_z^w$, $G_N^w$ and $G_E^h$ from the hydrostatic coefficients $L_z^h$, $G_N^h$ and $G_E^h$ in Eq. (1). But this is easier said than done, as the functions $m_h\left(e_0\right)$ and $m_w\left(e_0\right)$ in Eq. (1) have almost the same dependence on the elevation angle $e_0$ (see Figure 1).

Figure 1.

The mapping functions $m_h$ (blue) and $m_w$ (red) plotted against each other for a typical GNSS station in Beijing (latitude: 39.6086$°$ N, longitude: 115.8922$°$ E), winter time on January 16th, 2012, with the VMF1 model [13], parameterized by inputting data from the ECMWF numerical weather model EAR-40 [31].

2. Basic assumptions at the core of the definition of mapping functions and gradients

Mapping functions, as they were introduced by Marini [6, 32] are based on the assumption of a totally layered atmosphere. This means that the refractivity n is only a function of height (the exact meaning of the word height is related to the definition of geoid). The ray equation of radio waves (including light) obeys, in the spherical approximation and again for a totally layered atmosphere (dependence on geocentric radius r of the refractivity n), the prime integral relation

$$\begin{array}{c}n\left(r\right)rcos\left(\mathrm{e}\right)=n\left(r_0\right)r_{0}cos\left(e_0\right)\end{array}$$E4

where $r$ is the geocentric radius, $r_0$ is the geocentric radius at the receiver location, $n\left(r\right)$ is the refractivity at geocentric radius $r$, $e$ is the angle between the tangent to the bended ray and the local horizon (the plane perpendicular to the direction of $r$ at height $r$). $e_0$ is the elevation angle of the tangent of the bended ray at the receiver location.

The details of the computation of the ray path can be found in [6, 33, 34]. The refractivity of the atmosphere is a function of pressure, temperature and water vapor contents. A formula widely used is the Smith and Weintraub formula [35], derived for laboratory conditions (air perfectly mixed), as

$$\begin{array}{c}\left(n-1\right)=K_1\frac{P_d}{T}+K_2\frac{e}{T}+K_3\frac{e}{T^2}\end{array}$$E5

where $P_d$ is the partial pressure of dry air in millibars, T is the temperature in Kelvin, e is the partial pressure of water vapor. $K_1,K_2\text{and}{K}_3$ are constants. The $P_d$ term corresponds to the “dry” part of the refractivity, the $e$ terms correspond to the “wet” part of the refractivity. Many authors have improved the coefficients $K_1,K_2$ and $K_3$ year after year [15, 36, 37].

This formula can be easily rewritten as

$$\begin{array}{c}\left(n-1\right)=K_1^{\prime }\frac{P}{T}+K_2^{\prime }\frac{e}{T}+K_3\frac{e}{T^2}\end{array}$$E6

Where $P=P_d+e$. This rewriting, was the first term is denominated as the hydrostatic component of the refractivity, was proposed by Davis et al. [7] and then has been widely accepted, but lead to a track of confusion in the literature between the meaning of “hydrostatic” and “dry”. The word “hydrostatic” has specifically no meaning in Eq. (6), other than indicating that the total pressure is used instead of the partial pressure of the non-wet (dry) air, as in Eq. (5). The word “hydrostatic” has a precise meaning in numerical weather models [38], where it indicates that the equilibrium of an air column is a balance between the vertical pressure gradient and the buoyancy forces, neglecting convective processes [39] as a simplification of the Navier–Stokes primitive Equations [40]. This is also the assumption made in the Saastamoinen model of the atmosphere propagation delays [41], with the total pressure P at ground level taken as a parameter (and with also the assumption of an atmosphere “at rest”).

To a good degree of approximation, the refractivity of air obeys a twofold exponential formula [42].

$$\begin{array}{c}n\left(r\right)=1.+\delta n_h+\delta n_w=1.+N_{h}exp\left(\frac{r-r_0}{H_h}\right)+N_{w}exp\left(\frac{r-r_0}{H_w}\right)\end{array}$$E7

The terms $N_h$, $H_h$ and $N_w$, $H_w$ have, respectively, a value of 250 $·$ 10⁻⁶, 8.7 km, 128 $·$ 10⁻⁶ and 2.7 km for the location of our geodesy observatory in Tahiti (from the fit of radiosounding data over a typical year). The scale height $H_w$ varies from 1.5 km to up to 8 km from place to place and according to a seasonal cycle [43]. For all practical GNSS purposes, one can consider that the water vapor is concentrated in the troposphere (from 8 km over the poles to 18 km at the Equator [44, 45], and that the atmosphere extends up to 100 km [46, 47]. The International Union of Telecommunications [48] recommends the use, for radio-link purposes, on a worldwide basis and for altitudes taken from sea level, of the formula (7), with $N_h$ = 315 $·$ 10⁻⁶, $H_h$ = 7.35 km, the wet part being omitted (it is in fact included as a worldwide average in $N_h$ and $H_h$).

The prime integral (4) allows two things: 1/the computation of the path, 2/the computation of the time delay along the path as

$$\begin{array}{c}L={\int }_{\mathit{\text{path}}}\mathit{n ds}\end{array}$$E8

The extra delay (in equivalent length) caused by the atmosphere is

$$\begin{array}{c}\mathit{\delta L}={\int }_{\mathit{\text{path}}}\left(n-1\right)\mathit{ds}\end{array}$$E9

By inserting Eq. (6) into Eq. (9) we get the separation of $\mathit{\delta L}$ into additive “hydrostatic” ${\mathit{\delta L}}_h$ and “wet” ${\mathit{\delta L}}_w$ delays. The ratios of ${\mathit{\delta L}}_h$ and ${\mathit{\delta L}}_w$ with respect to the corresponding values taken along a vertical path are by definition (as in Eq. (1)) the hydrostatic ($m_h$) and wet ($m_w$) mapping functions that only depend on the elevation angle $e_0$ of the tangent of the bended ray at the receiver location.

Davis et al. [10] pushed the physical analysis of Eq. (9) a little bit further by introducing the notion of gradients. This notion is also based on the basic assumption of a main dependence of the refractivity with respect to height, with the refractivity in the neighborhood of the receiver written as

$$\begin{array}{c}n=n_V\left(r\right)+\mathit{\text{small lateral terms}}\end{array}$$E10

where r is taken along the local vertical of the receiver, and $n_V$ is the variation of n along the vertical of the observation site (the value of $n$ at the receiver station is $n\left(r_0\right)=n_V\left(r_0\right)$). One can note that this writing violates, on a pure mathematical ground the dependence of n on only the geocentric radius, that was assumed for the computation of the path in Eq. (4) (i.e. no small lateral terms should be present). If we define a local frame with units vector $\left(\widehat{x},\widehat{y}\right)$ in the tangent plane perpendicular to the vertical direction of the station (usually defined by the North and East directions as in Eq. (1), we get, with also the assumption of a “flat Earth”, the approximation

$$\begin{array}{c}n\left(r;x,y\right)\simeq n_V\left(r\right)+{\left(\frac{\partial n}{\partial x}\right)}_{\left(r\right)}x+{\left(\frac{\partial n}{\partial y}\right)}_{\left(r\right)}y\end{array}$$E11

This is nothing else than a Taylor series, meaning that $x$ and $y$ are assumed to be small, and the subscript $\left(r\right)$ emphasizes that the partial derivatives of $n$ are varying with the height $r$ (i.e. they are not taken at $r=r_0$). For low elevation angles of the path, $x$ and $y$ are by no means “small”, and can reach up to several hundreds of kilometers. We can define Eq. (11) as a “cylindrical” expansion of the refractivity.

If we insert this in Eq. (9), we get

$$\begin{array}{c}\mathit{\delta L}\left(e_0\right)={\int }_{\mathit{\text{path}}}\left(n_V-1\right)\mathit{ds}+{\int }_{\mathit{\text{path}}}{\left(\frac{\partial n}{\partial x}\right)}_{\left(r\right)}x\mathit{ds}+{\int }_{\mathit{\text{path}}}{\left(\frac{\partial n}{\partial y}\right)}_{\left(r\right)}y\mathit{ds}\end{array}$$E12

If we now divide the first right term of Eq. (12) by

$$\begin{array}{c}\mathit{\delta L}\left(e_0\right)={\int }_{\mathit{\text{vertical}}}\left(n_V-1\right)\mathit{ds}\end{array}$$E13

We get

$$\begin{array}{c}\mathit{\delta L}\left(e_0\right)=m\left(e_0\right){\int }_{\mathit{\text{vertical}}}\left(n_V-1\right)\mathit{ds}+{\int }_{\mathit{\text{path}}}{\left(\frac{\partial n}{\partial x}\right)}_{\left(r\right)}x\mathit{ds}+{\int }_{\mathit{\text{path}}}{\left(\frac{\partial n}{\partial y}\right)}_{\left(r\right)}y\mathit{ds}\end{array}$$E14

where $m\left(e_0\right)\approx \frac{1}{sin{e}_0}$ is by definition the mapping function. The value $\frac{1}{sin{e}_0}$ is obtained by setting all the coefficients $a,b,c\dots$ to $0$ in Eq. (2).

By writing $R^2=x^2+y^2$, $x=Rcos\varphi$, $y=Rsin\varphi$, and taking advantage of the fact that the path is nearly a straight line, as $n$ is close to 1 at a 10⁻³ level, we can write, for the two integrals involving the derivatives of $n$, $R=\mathit{\text{rcotg}}\left(e_0\right)$ and $\mathit{ds}=\frac{\mathit{dr}}{sin\left(e_0\right)}$. This is permissible, because physically these derivatives, as well as x and y are assumed to be small quantities. We obtain for the integral relative to the partial derivative $\left(\frac{\partial n}{\partial x}\right)$

$$\begin{array}{c}{\int }_{\mathit{\text{path}}}{\left(\frac{\partial n}{\partial x}\right)}_{\left(r\right)}x\mathit{ds}=m\left(e_0\right)cot{e}_{0}cos\varphi {\int }_{r=r_0}^{r_{\mathit{top}}}{\left(\frac{\partial n}{\partial x}\right)}_{\left(r\right)}r\mathit{dr}\end{array}$$E15

where $r_{\mathit{top}}$ is the top of the atmosphere with respect to the geocentric radius (around 100 km), and a similar expression in $sin\varphi$ for the partial derivative $\left(\frac{\partial n}{\partial y}\right)$.

The precise details of the mathematical machinery linking Eq. (11) to Eq. (1) can be found in Davis et al. [10]. The important fact, from a physical point-of-view is that, if we split the refractivity into a “hydrostatic” and a “wet” part, we get the “hydrostatic” and “wet” gradients of Eq. (1) as

$$\begin{array}{c}{G}_N^h={\int }_{r=r_0}^{r_{\mathit{top}}}{\left(\frac{\partial n_h}{\partial x}\right)}_{\left(r\right)}r\mathrm{dr},G_E^h={\int }_{r=r_0}^{r_{\mathit{top}}}{\left(\frac{\partial n_h}{\partial y}\right)}_{\left(r\right)}r\mathrm{dr}\end{array}$$E16

$$\begin{array}{c}{G}_N^w={\int }_{r=r_0}^{r_{\mathit{top}}}{\left(\frac{\partial n_w}{\partial x}\right)}_{\left(r\right)}r\mathrm{dr},G_E^w={\int }_{r=r_0}^{r_{\mathit{top}}}{\left(\frac{\partial n_w}{\partial y}\right)}_{\left(r\right)}r\mathrm{dr}\end{array}$$E17

The significations of the gradients are therefore the integration, along the altitude, weighted by the altitude, of the North and East directional derivatives of the “hydrostatic” and “wet” parts of the refractivity, evaluated along the vertical of the receiver location. It is in fact an integration along the geometrical line-of-sight.

3. Physical meaning of zenithal delays and gradients

The modeling of the extra-delays caused by the atmosphere by the combination of mapping functions and gradients of Eq. (1) has proved very effective since Davis introduced his formula 30 years ago [49, 50, 51]. But what is the real meaning of effective?

We have to remember that this model was primarily introduced to model atmospheric delays in VLBI, then to improve positioning estimates from GNSS data, and it is now battle-proven for these two applications. But another application, being known today as GNSS meteorology, emerged during the nineties, first with the modeling of the integrated water vapor contents along the vertical of the GNSS receiver (i.e. no gradients), known as “precipitable water” (or PW), that used the $L_z^w$ zenithal delay converted to PW through a multiplicative constant, known as the $\mathrm{\Pi }$ constant introduced by Bevis et al. [52]. Because the wet and dry mapping functions cannot be separated, for any practical purposes, in Eq. (1), the separation between the sum $L_z^d+L_z^w$ and $L_z^w$ must be done by introducing an “external hydrostatic estimate” $L_z^h$, the model of choice being the so-called Saastamoinen model [41]. By its own inception, a PW time series is relative to a particular GNSS station, and does not provide any information about the lateral gradients of the water vapor contents of the atmosphere for this site. But a dense network of GNSS receivers do. An even more powerful way to grasp the 3D and even 4D (with the inclusion of time) variations of the water vapor contents of the atmosphere is the tomography, first promoted by [1, 53, 54]. In the approach of tomography, Eq. (1) is just seen as an intermediate tool, the data inputted in the tomography software being the reconstructed $\delta L^w$ (the “wet” part of Eq. (1)). The tomography approach needs a dense network of GNSS receivers over a limited area, and take advantage of a multiple crossing paths between the receivers and the satellites of the GNSS constellations to invert the intrinsically ill-posed correspondence between the $\delta L^w$ and the 3D atmospheric water vapor refractivity field over the area. All the tomography software treat, to obtain a tractable problem, the rays as straight lines. This means that low-elevation slant delays cannot be considered.

Some authors [51, 55, 56] tried to assess the physical meaning of tropospheric gradients, but their effort were limited to qualitative assessments and correlations studies. Up to our knowledge [57], nobody is using gradients as data to constraint operational NWMs, albeit efforts having made to extract gradients from NWM numerical simulations [14] or make comparisons with NWMs outputs [58], or even to propose the use of slant delays for such a use [59]. The only GNSS data products that are currently inputted (assimilated) in NWMs are total zenithal delays (i.e. the sum $L_z^d+L_z^w$), as for example in the latest Météo-France AROME model [60].

This is clearly sending the message that the meteorology community does not yet consider gradients as a usable data set. We think that the main reason for this is the underlying assumption of the cylindrical Taylor’s expansion [Eq. (11)], at the basis of the notion of gradients, where a strict separation between vertical variations and lateral variations is assumed, and supposed valid over all the troposphere (at least as seen from the receiver location). This assumption is closely related to the hydrostatic assumption, itself closely linked to the highly non-linear Navier–Stokes equations, which admit as solutions a combination of laminar and turbulent/convective flows. At scales larger than a few tens of kilometers, the atmospheric flows are mostly horizontal [61]. This corresponds to the highest resolution available for typical MNW models, built around the hydrostatic assumption [62]. The atmospheric turbulence [63] itself is organized as “vortices”, or eddies, with scales varying over several orders of magnitude, from a few meters to several hundreds of kilometers [64, 65]. A combination of laminar and turbulence is also possible, and it is known as “frozen flow”, where “frozen turbulence” is carried by laminar flow [66]. This is illustrated for the layman by clouds driven by the wind. Atmospheric turbulence/convection is modeled through statistical tools, the structure functions [67], that obeys an exponential decay with altitude (i.e. turbulence is “higher” in the boundary layer) [68]. The definition of gradients by Davis et al. [10] is simply too crude from a “meteorological” point-of-view.

4. Beyond zenithal delays and gradients

Therefore, what can be the future of the modeling of neutral delays in GNSS meteorology? Applications in GNSS positioning and VLBI clearly show that Eq. (1) is sufficient for these applications, because what is of interest to these users are the integrated delays, not directly the variations of refractivity in the atmosphere. Eq. (1) is sufficient by itself to model these slant (extra) delays, as evidenced by tomography applications and the statistical analysis of these delays [69]. The zenithal total delays have proven to have a physical meaning, as they are related to the modeling of PW through an a priori model of the “dry” atmosphere and a proportional correspondence to zenithal wet delays. They are also feeding the current medium resolution NWM models. The gradients themselves are more questionable. They are merely ad’hoc, empirical corrections introduced for positioning and VLBI applications.

Can the definition of gradients be improved? From a physical point-of-view, we do not think so. The main assumption to derive the delay gradients in Davis et al. formula (Eqs. (16) and (17)) is an integration, along the line-of-sight receiver-satellite, of the gradients of the refractivity. Even with a better “geometrical definition” of the gradients, taking into account the curvature of Earth, the bending of the rays, etc.…, the main problem is that a line-of-sight station-satellite usually cross – and average– many eddies. According to [70], the shape and size of the eddies depend on the altitude. Close to the ground (0–2 km of altitude), the eddies are assumed to be small and not far from isotropic, while the irregularities at higher altitudes are highly anisotropic, i.e., the eddies become more flattened laterally. Along the vertical, the refractivity variation is mainly dominated by an exponential decay [71], but this is not the case along the horizontal direction. Besides, the repartition of the lines-of-sight in the sky can be scarce or uneven. For example, the GPS constellation, the most used one because of the high quality of its orbit modeling, offer quasi-repeating repeating tracks where only a few satellites (4 to 12) are visible from a particular location (Figure 2). This means that only a few lines-of-sight can be used at any time, and that there is, from a practitioner point of view, not enough data to constraint a better representation of the slant delays than the six-parameters Eq. (1).

Figure 2.

The sky-tracks (in elevation and azimuth) of the GPS satellites (one color per satellite) visible from the THTI station (latitude: 17.5769$°$ S, longitude: 149.6063$°$ W), in the wet season on January 10th, 2018.

Hopefully, Augmented Constellations and Low-Earth-Orbits constellations (LEO) will become soon a reality [72, 73, 74], thanks to the ever-decreasing size and costs of satellites, as well as the availability of miniaturized atomic clocks [75]. LEO constellations are particularly interesting for GNSS meteorology, as their satellites will cross the sky in a few minutes instead of hours, with a boost by one order of magnitude, or even two, of the available line-of-sight geometries. Our proposal to keep the separation of the refractivity into a “hydrostatic” and “wet” part, with the “hydrostatic” slant part evaluated separately from proven models like the Saastamoinen [41] model and subtracted from the total slant delay, then to represent the wet refractivity field based on a mean exponential decay of the wet refractivity as

$$\begin{array}{c}\delta n_w\left(r\right)=N_{w}exp\left(\frac{r-r_0}{H_w}\right)\left(1.+{\epsilon }_w\left(x,y,z,t\right)\right)\end{array}$$E18

where the ${ϵ}_w$ terms represent the departure of the wet refractivity field from the exponential local decay law and $x,y,z,t$ are local coordinates with respect to a frame linked with the local GNSS receiver and $t$ is time. As the wet scale height can vary by a factor of four, it must be provided from external sources (for example from the ECMWF-ERA series of climate models, see [76]). An estimate of $H_w$ can also be determined from the slant wet delays themselves, but only if a reliable estimate of the wet refractivity is available, as the integral over the geometrical path between the GNSS satellite and the receiver is proportional to $N_w{H}_w$ for a pure exponential decay of the wet refractivity. Empirical relations also exist between the ground value of the refractivity and scale height for example [77], but they are probably out-of-date. $H_w$ is by itself a very important parameter, as [71] demonstrated that this scale height is related to the rate at which the PW decorrelates with horizontal separation.

On the contrary of Davis et al. [10], we fully represent the term ${\epsilon }_w\left(x,y,z,t\right)$ as a 3D (or 4D if the time is present) series expansion

$$\begin{array}{c}{\epsilon }_n\left(x,y,z,t\right)=\sum _n{\lambda }_n{\mathrm{\Phi }}_n\left(x,y,z,t\right)\end{array}$$E19

where the ${\mathrm{\Phi }}_n\left(x,y,z,t\right)$ are a set of suitably chosen orthogonal functions in the atmospheric lens comprised between the local horizon of the station and the local tropopause. The ${\lambda }_n$ are the coefficients of the expansion. If the shape of the tropopause boundary is known [78], the ${\mathrm{\Phi }}_n$ functions can be defined as empirical orthogonal functions (EOF) [79] or as a pre-defined set of orthogonal functions renormalized according to the Gram-Schmidt scheme [80].

A preliminary attempt with a small data set was made by [81] with the assumption of a constant altitude tropopause (see Figure 3), where the ${\mathrm{\Phi }}_n$ orthogonal functions are a subset of Zernike functions [82]. The line-of-sight are assumed to be straight-lines to obtain tractable equations, as it is the case for tomography [83, 84] and the statistical analysis of the slant delays [85, 86]. This implies that low-elevation rays cannot be taken into account.

Figure 3.

The geometry of the inversion of the wet delays, with the representation of eddies in the troposphere, flattened with altitude and pushed by the wind [62, 63].

The integral relation to be solved with respect to ${ϵ}_w$ is therefore

$$\begin{array}{c}\delta L_w\left(e_0\right)=N_w{\int }_{\begin{array}{c}\mathit{\text{geometrical}}\\ \mathit{\text{path}}\end{array}}exp\left(\frac{r-r_0}{H_w}\right)\left(1.+{\epsilon }_w\left(x,y,z,t\right)\right)\mathit{ds}\end{array}$$E20

This integral relationship is averaging the wet refractivity field along the lines-of-sight (fan-beam tomography [87, 88]), and the inversion in terms of ${\lambda }_n$ coefficients must be regularized. By construction, the ${\epsilon }_w$ correction must be small, so we can use a truncated Singular Value decomposition (the EOF approach) or a Tikhonov approach [89] to enforce this smallness with respect to 1. The use of a priori refractivity values along the vertical for sites collocated with radiosoundings can also be envisaged [90] (in preparation). The Tikhonov approach, and its ability to model local variations of the refractivity field has been investigated in the framework of radar tomography [87, 91, 92].

The only case where the hypothesis of a small ${ϵ}_w$ can be violated occurs during inversion episodes, where atmospheric temperature increases when altitude increases. The warm inversion layer then acts as a cap and stops atmospheric mixing [93] causing a large deviation of the refractivity with respect to the exponential decay.

The end-product for the meteorology community of the inversion of Eq. (19) cannot only be the set of $H_w$ and ${\lambda }_n$ coefficients, that are too difficult to handle. We propose, in addition, to give the results in the form of records over a grid the resolution of which is in agreement with the maximum degree of the expansion in Eq. (19), with respect to a suitable ellipsoid (like WGS84), and with these fields:

Observation Time, Latitude, Longitude, Geometrical height, total refractivity, wet refractivity.

The refractivity fields can then be converted, if needed, to water vapor levels according to Eq. (6) with suitable temperature profiles over the troposphere and/or feed high resolution NWM taking natively into account turbulent/convective processes [94]. Xia et al. [95] tried to derive the refractivity field from slant delays by substituting Eq. (6) into Eq. (9), but the underlying hypothesis is an atmosphere at rest, in a similar fashion of the neutral delay model of Saastamoinen [41, 96].

Is the approach developed in this article directly implementable in GNSS software, as a replacement of the usual approach of Eq. (1)? The response is a careful yes [69]. Strictly speaking, a mapping function defines, from the point of view of differential geometry, a time-evolving coordinate chart that is a non-orthogonal system of coordinates made of the refracted elevation at ground level, the length along the bended ray, and the azimuth. We think that such an implementation in GNSS software implies at least the use of a constant (i.e., not evolving with time) system of coordinates (i.e., a constant mapping function), that therefore must be computed with respect to some standard model of the atmosphere, carefully designed and normalized [97]. For this purpose, it should be noted that the variation of the propagation delay caused by the bending is of second order with respect to the integration of the refractivity along the path [98].

Finally, the modeling of the wet refractivity field through an expansion series in time and space (Eq. (19)) can be also used to model tropospheric delays, in a correlated way, between uplink and downlink signals to planetary space crafts, where the uplink and downlink separation in time can reach tens of minutes or even hours [99].

5. Conclusion

We discussed in this brief paper the pros and cons of the standard approach mapping functions + gradients to model the neutral delays of the atmosphere, and more specifically the wet delays caused by the presence of water vapor in the troposphere. If this standard approach is almost perfect for people doing positioning, deformation and VLBI studies, as they see the neutral delays as “noise”, it is not so well adapted to people looking at these delays as signals to study atmospheric processes. In particular, the standard definition of gradients is too crude, and does not permit to have access to the horizontal turbulence/convection scales, that are key parameters to model these processes in high resolution NWM models. We therefore propose an alternative way to model the wet tropospheric delays, through a representation of the wet refractivity field as a perturbation over an exponential decay with altitude with a locally adjusted scale height and a time/space series expansion over a suitable basis of orthogonal functions. Our approach is computationally expensive, and maybe not suited for real-time applications, but its end-product are records of the total and wet refractivity values with high-resolution in time (minute-scale) and distance (sub km-scale), in accordance with the needs of future numerical weather models [38], the emerging field of the modeling of atmospheric rivers [100, 101] and besides does not require the additional step of water vapor tomography, with lower cost, better mobility and simpler operation [102].

Acknowledgments

This research was funded by a DAR grant from the French Space Agency (CNES) to the Geodesy Observatory of Tahiti. Feng Peng did this research in the framework of his Ph.D., with a 14-months stay at the Observatory of Tahiti funded through the Cai Yuanpei program.

Conflict of interest

The authors declare no conflict of interest.