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.
- GNSS meteorology
- deep space tracking
- neutral delays
- mapping functions
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 .
where is the slant (extra) delay with respect to propagation in vacuum along the bended ray, and, are the hydrostatic and wet zenith delays, and are the satellite elevation and azimuth angles as seen from the station, respectively; and , and are the north and east components of the hydrostatic and wet delays gradients; and 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 (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) and wet delay , essentially caused by the water vapor. The mapping functions are usually written  in the form of continuous fractions, that were introduced by Marini [6, 7] and normalized by Herring  in the form
Other (simplified) forms of the mapping functions can be found in the literature , but the mainstream form is always Eq. (2). The gradients themselves, noted by the upper case letter 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
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 ), with respect to the total zenithal delay (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  is therefore of the uttermost importance, as evidenced by the so-called Energy Balance models [17, 18] that can be written
Where is the solar constant (1360 W/m2),
The coefficients and are albedos, respectively in the visible and infrared wavelengths, both mainly driven by the water vapor contents of the Earth atmosphere . The coefficients reflects the cloud coverage, typically today at the 30% level, and the coefficient 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 means a cooling of Earth surface, and an increase of means a warming, with a lot of intricacies between the positive and negative feedbacks related to the water vapor cycle of the climate models . The ultimate goal of global long-term climate models  is to predict which effect will prevail (this is the
The study  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 , 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 , photometers , and water vapor radiometers . The only source providing
It is therefore important to separate the water vapor modeling coefficients , and from the hydrostatic coefficients , and in Eq. (1). But this is easier said than done, as the functions and in Eq. (1) have almost the same dependence on the elevation angle (see Figure 1).
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
where is the geocentric radius, is the geocentric radius at the receiver location, is the refractivity at geocentric radius , is the angle between the tangent to the bended ray and the local horizon (the plane perpendicular to the direction of at height ). 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 , derived for laboratory conditions (air perfectly mixed), as
where is the partial pressure of dry air in millibars,
This formula can be easily rewritten as
Where . This rewriting, was the first term is denominated as the hydrostatic component of the refractivity, was proposed by Davis et al.  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 , 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  as a simplification of the Navier–Stokes primitive Equations . This is also the assumption made in the Saastamoinen model of the atmosphere propagation delays , with the total pressure
To a good degree of approximation, the refractivity of air obeys a twofold exponential formula .
The terms , and , have, respectively, a value of 250 10−6, 8.7 km, 128 10−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 varies from 1.5 km to up to 8 km from place to place and according to a seasonal cycle . 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  recommends the use, for radio-link purposes, on a worldwide basis and for altitudes taken from sea level, of the formula (7), with = 315 10−6, = 7.35 km, the wet part being omitted (it is in fact included as a worldwide average in and ).
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 into additive “hydrostatic” and “wet” delays. The ratios of and with respect to the corresponding values taken along a vertical path are by definition (as in Eq. (1)) the hydrostatic () and wet () mapping functions that only depend on the elevation angle of the tangent of the bended ray at the receiver location.
Davis et al.  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
This is nothing else than a Taylor series, meaning that and are assumed to be small, and the subscript emphasizes that the partial derivatives of are varying with the height (i.e. they are not taken at ). For low elevation angles of the path, and 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
If we now divide the first right term of Eq. (12) by
where is by definition the mapping function. The value is obtained by setting all the coefficients to in Eq. (2).
By writing , , , and taking advantage of the fact that the path is nearly a straight line, as is close to 1 at a 10−3 level, we can write, for the two integrals involving the derivatives of , and . This is permissible, because physically these derivatives, as well as
where is the top of the atmosphere with respect to the geocentric radius (around 100 km), and a similar expression in for the partial derivative .
The precise details of the mathematical machinery linking Eq. (11) to Eq. (1) can be found in Davis et al. . 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 zenithal delay converted to PW through a multiplicative constant, known as the constant introduced by Bevis et al. . Because the wet and dry mapping functions cannot be separated, for any practical purposes, in Eq. (1), the separation between the sum and must be done by introducing an “external hydrostatic estimate” , the model of choice being the so-called Saastamoinen model . 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 (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 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 , nobody is using gradients as data to constraint operational NWMs, albeit efforts having made to extract gradients from NWM numerical simulations  or make comparisons with NWMs outputs , or even to propose the use of slant delays for such a use . The only GNSS data products that are currently inputted (assimilated) in NWMs are total zenithal delays (i.e. the sum ), as for example in the latest Météo-France AROME model .
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 . This corresponds to the highest resolution available for typical MNW models, built around the hydrostatic assumption . The atmospheric turbulence  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 . This is illustrated for the layman by clouds driven by the wind. Atmospheric turbulence/convection is modeled through statistical tools, the structure functions , that obeys an exponential decay with altitude (i.e. turbulence is “higher” in the boundary layer) . The definition of gradients by Davis et al.  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 . 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
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 , 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 , 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 . 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  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 terms represent the departure of the wet refractivity field from the exponential local decay law and are local coordinates with respect to a frame linked with the local GNSS receiver and 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 ). An estimate of 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 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 , but they are probably out-of-date. is by itself a very important parameter, as  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. , we fully represent the term as a 3D (or 4D if the time is present) series expansion
where the 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 are the coefficients of the expansion. If the shape of the tropopause boundary is known , the functions can be defined as empirical orthogonal functions (EOF)  or as a pre-defined set of orthogonal functions renormalized according to the Gram-Schmidt scheme .
A preliminary attempt with a small data set was made by  with the assumption of a constant altitude tropopause (see Figure 3), where the orthogonal functions are a subset of Zernike functions . 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.
The integral relation to be solved with respect to is therefore
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 coefficients must be regularized. By construction, the correction must be small, so we can use a truncated Singular Value decomposition (the EOF approach) or a Tikhonov approach  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  (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 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  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 and 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:
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 . Xia et al.  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 . 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 . 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 .
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 .
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 , 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 .
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.