On horizontal distribution of vertical gradient of atmospheric refractivity

The horizontal distribution of the vertical gradient of atmospheric refractivity is important for the assessment of the propagation of electromagnetic waves on nearly horizontal paths. A 5‐year set of meteorological data, obtained from the ERA‐Interim numerical weather product, has been used to analyze this distribution statistically. Vertical gradient maps of the Central European region have been processed to derive empirical probability distributions of the difference between local point gradient values at two locations. The difference of the effective (path‐averaged) gradient along the path and the local gradient is shown to be statistically distributed so that the quantiles increase linearly with distance in the interval from 100 to 1000 km. Finally, a spatial correlation function is obtained and described by an exponential model with correlation distances determined in the range of 400–700 km.


Introduction
The atmosphere, as a dielectric medium, can be characterized by the spatial and temporal distribution of the refractive index of air that is dependent on pressure, temperature and humidity. Since the atmosphere is predominantly horizontally stratified so is the refractive index and it is the vertical component of the gradient which is the most important parameter describing the inhomogeneity of the atmosphere as a dielectric medium.
Spatial distribution of atmospheric refractivity is known to affect the propagation of electromagnetic (EM) waves in the atmosphere (Kerr, 1987). In particular, propagation on nearly horizontal paths is directly affected by the vertical gradient of refractivity due to the effect of the bending of EM waves. Statistics of the gradient obtained from long-term atmospheric measurements (e.g. Hall and Comer, 1969;Akiyama, 1977;Grabner and Kvicera, 2006) are used to predict statistical propagation parameters (ITU, 2015).
In frequency bands above 1 GHz that are widely utilized in terrestrial communication systems, typical path lengths are usually limited to less than 100 km. In this case, the vertical gradient of refractivity is often assumed to be constant along the whole propagation path and single-point statistics of the gradient are applied to assess the long-term link performance (ITU, 2015). On the other hand, there are applications such as unmanned aerial vehicle systems or low elevation links (Vanhoenacker-Janvier et al., 2013) with transmitters in higher altitudes operating on significantly longer paths. Here the knowledge about the horizontal variation of the vertical gradient is needed in order to predict propagation on such paths accurately. An appropriate statistical modeling of horizontal spatial distribution of the vertical gradient can be applied in remote sensing too (Jicha et al., 2013).
In this study, the two-dimensional spatial distribution of the vertical gradient of atmospheric refractivity is analyzed based on 5-year meteorological data. The horizontal variation of the gradient is described by statistics of gradient difference and path-average gradient along the path. Finally, the spatial correlation characteristics of the gradient are derived and modeled.

Input data and processing
Numerical weather products ERA-40 and, recently, ERA-Interim (ERAI) (Dee et al., 2011) by the European Centre for Medium-Range Weather Forecasts (ECMWF) are regularly used in ITU-R SG 9 as a source of long-term and global meteorological data for the propagation modeling purposes. The advantages of the ERAI dataset are twofold: (1) ERAI dataset is carefully checked by ECMWF and is free of serious errors and (2) it provides global coverage with sufficient spatial resolution 0.75 ∘ and temporal resolution 6 h. The ERAI data has been processed recently in order to develop global statistical maps of atmospheric refractivity parameters (Grabner et al., 2014). Spatial correlation characteristics of the vertical refractivity gradient have been initially studied in (Grabner et al., 2015). The vertical profiles of atmospheric refractivity were extracted from the ERAI database of meteorological parameters and vertical gradients were calculated by linear regression.
The following parameters were extracted from the ECMWF ERAI database with a time resolution of 6 h: Horizontal distribution of vertical gradient of atmospheric refractivity 295 model level temperature T (K) in 12 model levels corresponding to the lowest part of the atmosphere (levels 49-60), model level specific humidity q (kg kg −1 ) in the same 12 model levels, and surface pressure p s (Pa). Additional time invariant parameters were also obtained, such as surface geopotential Φ s (m 2 s −1 ).
The pressure values p i at model level boundaries (i = 1, … , 61) are derived from the surface pressure p s . The heights h i above the surface corresponding to a particular model levels with the pressure p i are then obtained from the geopotential values Φ i . Those are, in turn, derived from Φ s using a recommended standard iterative procedure that is determined from the discrete analog of the hydrostatic differential equation (ECMWF, 2016).
The atmospheric refractivity N (N-unit) (N = (n − 1)10 6 , n is the refractive index of air) for radio waves (0-300 GHz at least) is dependent on pressure p (hPa), temperature T (K) and partial water pressure e (hPa). The relationship is often expressed by the formula: where for the Rüeger best-average model (Rüeger, 2002), K 1 = 77.689 ± 0.0094, K 2 = 71.295 ± 1.3, K 3 = (3.75463 ± 0.0076)10 5 . Partial water pressure e is obtained from the specific humidity of air q by the relation: where R dry = 287.0597 J kg −1 K −1 and R vap = 461.5250 J kg −1 K −1 are the gas constants of dry air and water vapor. The vertical gradient of refractivity G (N-units km −1 ) is derived by means of the linear regression of refractivity values at heights below the specified limit h max . Thus, the vertical profiles of refractivity values N i (h i ) for heights h i < h max were approximated by a linear equation: where h denotes the height above the surface and N (0) is an additive constant. In the following, h max = 100 m above the ground is chosen because, while refractivity at altitudes higher than 100 m is also relevant, the focus on the lowest troposphere is motivated by the higher spatial variability of refractivity near the ground and by propagation being more affected by the refractivity of the lowest layer. The heights obtained from the model levels as described above are not equidistant. For example, in Prague they are about: 0, 23, 58, 107, 173, 258, … m above the ground. In this case, the refractivity at 100 m is approximated first by linear interpolation (from heights 58 and 107 m) and then the gradient is obtained by a linear regression of refractivity at heights 0, 23, 58 and 100 m above the ground. An example of the spatial distribution of the vertical gradient less than 100 m above the surface is shown in Figure 1. The spatial distributions depicted correspond to the times 0000 and 0600 UTC. The standard value of the gradient G 100m ≈ − 40 N -units km − 1 is present over most of continental Europe and North America but super-refractive gradients are present over some sea and coastal areas.

Point gradient
In order to demonstrate the spatial variability of the vertical gradient of refractivity, the gradient difference dG = G(0) − G(d) between two locations separated by a distance d is analyzed. Figure 2 shows the cumulative distributions (CDs) of dG for distances from 0 to 1000 km. The empirical distributions were calculated using data from the years 2008-2012 considering paths with a central location in Prague. The other points are given by eight azimuthal angles (measured from the northerly direction) ranging from 0 to 315 ∘ and by distances from 0 to 1000 km. The maximal distance (1000 km) analyzed was chosen as a practical upper bound since: (1) correlation is low on greater distances, see below, and (2) in many propagation scenarios, the propagation path of wanted or unwanted signals is likely shorter than this limit.
Note that Figure 2 shows both CDs and complementary is the usual CD) in order to clearly see both tails of the distribution. The difference dG, CD is slightly asymmetrical toward the negative differences. The gradient difference extremes are clearly increasing with the distance of the two locations.

Effective gradient
Considering propagation paths longer than approximately 20-50 km, point gradient is no longer the most suitable parameter for propagation estimation. Instead, the averaged propagation effect is better characterized by the effective gradient (Mojoli, 1980;ITU, 1996). The effective (or path-average) gradient is obtained as a mean value of the gradient along the propagation path: The effective gradient is usually more suitable to describe propagation on longer paths than the point gradient because local ray curvature is related to the local gradient and overall ray bending along the path is approximately given by a summation of local ray bending along the whole path.
From the propagation point of view, the extremal gradients occurring only rarely have the most significant impact on the propagation and their incidence statistics are of practical importance for the effective design of radio links. In this respect, time percentages as low as 0.01% are relevant. The effective gradient exceeded for 0.01% of time is obtained from ERAI data and compared with the statistics provided by ITU-R P.530 (ITU, 2015). Figure 2 of ITU-R P.530 shows the dependence of the effective k-factor k e exceeded 99.99% of the time on the path length for the continental temperate climate. Applying the known relation between the k-factor and the vertical gradient, the percentiles k e99.99 % are equivalent to the effective gradient percentiles according to the relation: ITU-R P.530 gives G e0.01 % (100) ≈ 20 N -units km − 1 for 100 km paths and G e0.01 % (200) ≈ 0 N -units km − 1 for 200 km paths. On the other hand, the 0.01% percentiles G e0.01 % obtained from ERAI for the location of Prague decrease from about 30 to 0 N-units km −1 for distances ranging from 100 to 1000 km. Therefore, for path lengths above 100 km, ITU-R P.530 and ERAI agree sufficiently well.
There is a higher inconsistency between ITU-R P.530 and ERAI statistics for path lengths less than approximately 100 km. But in this region, the accuracy of the ERAI results should be naturally limited since the resolution 0.75 ∘ is not sufficient to analyze accurately the length scales below 80 km. (One must also take into account that the ITU-R P.530 statistics of k e are originally derived indirectly from propagation measurements of a received signal level but the ERAI results come directly from meteorological data.) On the other hand, it is confirmed that G e quantiles are approaching a constant value with further increasing path lengths.
Clearer path length dependence is revealed when analyzing the difference between the point gradient and the effective gradient. Figure 3 shows the CDs of the difference dG e = G(0) − G e (d) for different path lengths d.
Both CDs and complementary CDs are depicted in Figure 3. The CD of the difference dG is highly symmetrical and with lower extremes than the point difference above which seems to be a logical consequence of the integral definition of G e . The difference dG e extremal values are clearly increasing with an increasing path length. For example, the 0.1 % quantiles are from about −40 to −125 N-units km −1 for path lengths from 100 to 1000 km. While the results presented in this section could find their main application in EM wave propagation problems, the horizontal distribution of the vertical gradient has not been extensively studied so far. An important first step of the gradient structure modeling is presented in Section 4 by a derivation of a spatial correlation of the spatial random field (SRF).

Spatial random field characteristics of a gradient
Assuming the spatial distribution of the refractivity gradient as SRF (Christakos, 2005), its statistical parameters were extracted from the analyzed dataset. The continuous SRF is generally described by a joint probability density function (PDF) of the dimension m for any m, but for practical reasons only the case m = 2 is usually considered when the empirical PDF is being extracted from the data. Figure 4 shows the examples of two-dimensional joint PDFs of the vertical gradients at locations separated by a surface distance d from 50 to 950 km with the central location in Prague. Clearly, the degree of correlation between gradient values is decreasing with the horizontal distance. The single, as well as joint, PDFs of the gradient are not generally of the Gaussian type and are difficult to model analytically.
[The single PDF of the gradient can be, e.g. expressed by the linear combination of three Gaussian distributions (Grabner and Kvicera, 2011).] Therefore, it is more convenient to use a correlation theory approach to model a horizontal structure of the gradient. The spatial correlation function (SCF) is defined as: where x 1 , x 2 denote different spatial locations separated by distance d = |x 1 − x 2 | and covariance C is obtained from: where G(x) is the gradient at location x and E{} denotes a mean value that is approximated by time averaging over the analyzed time period. Figure 5 shows the SCFs obtained for the central locations of Prague, London, New York and Denver (i.e. x 1 in Equations (6) and (7)) from the dataset (2008)(2009)(2010)(2011)(2012). The other locations x 2 were considered to be at distances d (km) from the central location with an azimuth angle from 0 to 315 ∘ . The correlation is seen as being fairly isotropic in Prague and Denver, but some anisotropy is evident in New York. Such anisotropy could be caused by the proximity of the ocean-land transition. However, the results from London show that it is not straightforward to relate the anisotropy to this effect in general. The mean correlation function averaged over azimuth angles, also depicted in Figure 5, decreases with distance in a way that can be described by the following exponential model: with the correlation distances d c = 631, 424, 430, 674 km for all four locations, respectively. These correlation distances were obtained by a least-square fitting of the model Equation (8) to the average (averaging over azimuths) correlation functions as shown in Figure 5. The root mean square error of the fitted model values f i with respect to the obtained average correlations is equal to 0.033, 0.038, 0.030 and 0.058, respectively. Note that the correlation distance d c is expected to be a site-specific parameter dependent on local climatic and orographic conditions. The SCF obtained is the simplest way to model spatial random distribution of the gradient. Using SCF allows, in principle, to model a joint PDF of a gradient or generate a SRF of the gradient from the SCF and from a single point PDF obtained from the local measurement.

Discussion and conclusions
The results presented show that on the horizontal distances within the interval (100, 1000 km), the difference between the local vertical gradients of refractivity at different locations is statistically distributed, so that its quantiles increase with the distance. The effective gradient along the propagation path that better characterizes propagation on extended paths is shown to be different from the local gradient in such a manner that the quantiles of the difference dG e increase with distance (or path length) almost linearly. These statistical results (CDs and joint PDFs) were determined for the climate region of Central Europe and they are therefore site-specific to some extent. Nevertheless, the qualitative conclusions presented were determined to be valid also for other inland locations. The SCF obtained is well-approximated by the exponential model Equation (8) with correlation distances between about 400 and 700 km depending on the location. It has been demonstrated that path direction may also influence the correlation distance at particular areas due to anisotropy of the gradient horizontal distribution.
Finally, propagation effects due to the horizontal distribution of refractivity gradient are not analyzed in this study. Interested readers can refer to the work of Valtr and Pechac (2005) who studied this problem using an analytical ray-tracing technique.