## 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.

### Keywords

- GNSS meteorology
- positioning
- VLBI
- deep space tracking
- neutral delays
- mapping functions
- gradients

## 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].

where

The mapping functions “map” the so-called slant neutral atmosphere (extra) delay

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

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

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

Where ^{2}), * T*is the mean temperature on the Earth surface in Kelvin,

*is the time.*t

The coefficients * 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

## 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

*of the refractivity*r

*), the prime integral relation*n

where

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

where * T*is the temperature in Kelvin,

*is the partial pressure of water vapor.*e

This formula can be easily rewritten as

Where * 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].

The terms ^{−6}, 8.7 km, 128 ^{−6} 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 ^{−6},

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

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

By inserting Eq. (6) into Eq. (9) we get the separation of

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

where * r*is taken along the local vertical of the receiver, and

*along the vertical of the observation site (the value 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*n

This is nothing else than a Taylor series, meaning that

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

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

We get

where

By writing ^{−3} level, we can write, for the two integrals involving the derivatives of * x*and

*are assumed to be small quantities. We obtain for the integral relative to the partial derivative*y

where

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

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

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

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).

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

where the

On the contrary of Davis et al. [10], we fully represent the term

where the

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

The integral relation to be solved with respect to

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

The only case where the hypothesis of a small

The end-product for the meteorology community of the inversion of Eq. (19) cannot only be the set of

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.