1.

Introduction

The effect of atmospheric refraction on light has been known and studied since antiquity.¹ Considerable work has been accomplished on refractive effects related to astronomical observations.^2–4 For terrestrial propagation scenarios, refractive effects can range from relatively subtle effects, such as an apparent shift in object position, to things as spectacular as mirages and the green flash.^5,6 Extreme events such as mirages and green flash have been studied extensively and their physics is well known.⁷ However, very few efforts have been directed at measuring and quantifying the amount of refractive bending one might expect on an ordinary day. With the development of the Global Positioning System (GPS) satellite network, refractivity profiles have been developed using radio occultation of GPS signals.⁸ In 2013, a patent was awarded to the Boeing Company for a method to measure the refractivity profile of a parcel of atmosphere using airborne imaging sensors.⁹

An imaging approach to measure the variations in the path-averaged refractive index gradient is described here. The sea level vertical refractive index gradient in the U.S. Standard Atmosphere model is $-2.7\times {10}^{-8}{\mathrm{m}}^{-1}$ at a wavelength of 500 nm. At any location, the actual gradient varies due to turbulence, local weather conditions, and time of day. It is important to understand the nature and extent of these variations and how they correlate to changing weather conditions. The imaging approach can be applied to any near-horizontal path through the atmosphere and provides a low-cost solution for ground truth measurement of refractive bending. From the refractive bending measurements, estimates of the refractive index gradient can be obtained. Observations can be done over several days, or months at a time, with the camera in time-lapse mode. Besides being low-cost, the ability to measure refractive bending at visible or infrared wavelengths is an added advantage. The imaging-based quantification of refractive bending can help to validate and improve fine-scale numerical weather prediction models by providing quantification of the vertical temperature lapse rate in the lower atmosphere. The vertical temperature lapse rate in the lower atmosphere is strongly influenced by the radiative properties of the surface, thus it can change significantly with cloudiness, terrain, and changes in ground cover (e.g., snow) yet it is only directly measured approximately twice a day with radiosonde balloon launches that are widely spaced geographically. It is currently difficult to capture the radiative nuances of partial cloudiness, terrain, and ground cover with numerical weather models. Thus if refractive bending observations were conducted with low-cost digital photography along with standard surface weather observations, the deduced temperature lapse rate could be used to increase forecast accuracy in the lower atmosphere.

An experiment was conducted at the Air Force Institute of Technology (AFIT) to demonstrate the time-lapse imaging approach. The methodology is described in Sec. 2. Section 3 discusses the results from the experiment and the correlation to meteorological conditions. A summary of the work and future research directions are presented in Sec. 4.

2.

Methodology

2.1.

Relationship Between Refractive Index Gradient and Apparent Change in Object Position

Figure 1 shows an example ray trajectory through the atmosphere. Due to refractive bending of the ray, to a remote observer, the object appears to be displaced from its actual position. From Snell’s law, it follows that, for near-horizontal paths through the atmosphere,¹⁰

Eq. (1)

$$\frac{\mathrm{d}n}{\mathrm{d}z}=-n\kappa ,$$

where $n$ is the refractive index of air, $\mathrm{d}n/\mathrm{d}z$ is the vertical refractive index gradient, and $\kappa$ is the curvature of the ray trajectory at that location.

Fig. 1

A ray trajectory through layered atmosphere and the associated shift due to refractive bending.

The refractive index is related to the air density, $\rho$, via the following equation:

Eq. (2)

$$n-1=a\rho ,$$

where $a$ is a wavelength-specific constant derived from the Edlen66 relation¹¹ for the dry air index of refraction normalized for an air density of $1\mathrm{kg}{\mathrm{m}}^{-3}$. Edlen’s equation,

Eq. (3)

$$n-1={10}^{-6}[83.4212+24060.3/(130-1/{\lambda }^2)+159.97/(38.9-1/{\lambda }^2)],$$

where $\lambda$ is the wavelength in $\mu \mathrm{m}$, actually calculates the refractivity (the portion within the brackets) that is dependent on wavelength and applicable in the UV to longwave infrared ($\sim 20\mu \mathrm{m}$). For visible light at 500 nm, $a$ in Eq. (2) is the $n-1$ value divided by $1.225\mathrm{kg}{\mathrm{m}}^{-3}$ (sea level air density in a standard atmosphere) and is equal to $0.227\times {10}^{-3}{\mathrm{m}}^3{\mathrm{kg}}^{-1}$. From the ideal gas law,

Eq. (4)

$$\rho =\frac{p}{R_{\mathrm{d}}{T}_{\mathrm{v}}},$$

where $p$ is the pressure, $T_{\mathrm{v}}$ is the virtual temperature, and $R_{\mathrm{d}}=287\mathrm{J}{\mathrm{kg}}^{-1}{\mathrm{K}}^{-1}$ is the dry air gas constant. The virtual temperature is a convenient way to quantify the effect of water vapor on air density by raising the dry air temperature, a proportional amount to represent the quantity of lower density water vapor contained in the ambient air. This is shown in the following equation:

Eq. (5)

$$T_{\mathrm{v}}=(1+0.61q)T,$$

where $q$ is the specific humidity (the ratio of the mass of the water vapor over the mass of the total air) and $T$ is the actual air temperature.

By using Eqs. (1), (2), and (4) and $\partial p/\partial z=-\rho g$, where $g$ is the acceleration due to gravity, a relationship between the curvature of an atmospheric ray path and vertical gradients of air density and temperature can be derived:

Eq. (6)

$$\kappa =-\frac{a}{n}\frac{\mathrm{d}\rho }{\mathrm{d}z}=\frac{n-1}{n{T}_{\mathrm{v}}}(\frac{g}{R_{\mathrm{d}}}+\frac{\partial T_{\mathrm{v}}}{\partial z}).$$

Thus the curvature can be directly calculated from the index of refraction and the vertical virtual temperature gradient. Typical values of $\partial T_{\mathrm{v}}/\partial z$ in daytime neutrally buoyant conditions are near $-9.8\mathrm{K}{\mathrm{km}}^{-1}$, but can be greater than zero in nighttime stable, temperature inversion conditions. Equation (6) dictates that typical nighttime conditions with a temperature inversion near the ground $(\partial T_{\mathrm{v}}/\partial z>0)$ will create more ray path curvature and objects near the horizon will appear higher at night than they will during the day when the lower atmosphere has a near neutrally buoyant lapse rate.¹²

Referring back to Fig. 1, if the slope of the ray is close to zero, which is true for near-horizontal paths,

Eq. (7)

$$\frac{{\mathrm{d}}^{2}z}{\mathrm{d}{x}^2}=-\kappa (x),$$

where $\kappa (x)$ is the curvature of the ray at any location $x$ along the path. By integrating Eq. (7) twice, the ray height at any location along the path can be determined:

Eq. (8)

$$z=-{\int }_0^x{\int }_0^{x_1}\kappa (x^{\prime })\mathrm{d}{x}^{\prime }\mathrm{d}{x}_1+(H-d)(x/L)+d,$$

where $L$ is the horizontal separation between the object and the observer, $H$ is the apparent shift in object position, and $d$ is the elevation of the observer with respect to the object. Since $z=0$ at $x=L$,

Eq. (9)

$$H={\int }_0^L{\int }_0^{x_1}\kappa (x^{\prime })\mathrm{d}{x}^{\prime }\mathrm{d}{x}_1={\int }_0^L(L-x)\kappa (x)\mathrm{d}x.$$

The second relationship in Eq. (9) was derived from the first by using integration by parts. Equation (9) shows that the curvature and, hence, the refractive index gradient is linearly weighted along the path, with maximum weight at the observer location and zero weight at the object location, i.e., variations in the index gradient at the object location do not affect the apparent object position. This is easily seen by recognizing that ray bending near the object will affect which ray an observer sees, but not the direction from which this ray arrives. Assuming $n\approx 1$, from Eq. (9), the average refractive index gradient along the path is

Eq. (10)

$${(\mathrm{d}n/\mathrm{d}z)}_{\mathrm{av}}=-{\kappa }_{\mathrm{av}}=-H/{\int }_0^L(L-x)\mathrm{d}x=-2H/L^2.$$

If the average refractive index gradient along the path changes by $\mathrm{\Delta }$, the apparent position of the object changes by $\mathrm{\Delta }{L}^2/2$. By measuring the position change, an estimate of the gradient change can be obtained.

2.3.

Stitching Image Shifts from Different Reference Images

Figure 2 shows example images taken during the afternoon and night of July 25, 2014. Due to different illumination conditions, the daytime images correlated very poorly with the night time images, thus making it difficult to estimate the relative shifts between daytime and nighttime images. Hence, two reference images were chosen: a midday image and a twilight image. The midday image correlated well with all the daytime images, but correlated poorly with the evening and nighttime images. The twilight image, on the other hand, showed good correlation with the late afternoon, evening, and nighttime images, but poor correlation with the morning and early afternoon images. In order to obtain the relative image motion, a cross-correlation algorithm was run between the reference images and every other image in the image bank. The relative shifts with respect to the two reference images are shown in Fig. 3. The jumps/discontinuities in the shift profiles suggest poor correlation. It is apparent from Fig. 3 that during certain periods, especially at daybreak and in the afternoons, the two profiles overlap well. While one profile starts getting noisy after the overlap zones, the other profile picks up right at these zones. The two profiles were stitched at the overlap areas, thus replacing the noisy information from one profile with the useful information obtained from the other profile.

Fig. 2

Example (a) daytime and (b) night time images. Note that the hospital appears to be more elevated in (b) than in (a).

Fig. 3

Relative image motion computed from the cross-correlation of all images with the midday and twilight reference images. Jumps/discontinuities represent poor correlation.

3.

Results

Figure 4 shows changes in the apparent hospital position from the evening of July 24, 2014, to the evening of August 3, 2014. The profile was obtained using the stitching technique described in Sec. 2.3. Gaps in the profile indicate cloudy/foggy/rainy periods when visibility was poor. Two components of image motion are evident in the profile. The faster, random motion is due to turbulence while the slow, vertical drift is due to changes in average refractive index gradient along the path. The thermal inversion that built up during the night caused the index gradient to decrease more with height and hence increase in the curvature of the rays coming from the hospital. During the day, as the ground and the surrounding air warmed up, the thermal inversion disappeared and the index gradient relaxed from the nighttime values and resulted in reduced ray curvature. The hospital building thus appeared to move up during the night and move down during the day. The image motion over the course of a day varied from about a meter on some days to about six times as much on certain days. Using Eq. (10), this corresponds to diurnal variations of $1.2\times {10}^{-8}$ to $7.2\times {10}^{-8}{\mathrm{m}}^{-1}$ in the average refractive index gradient. The gradient changes were more pronounced during the mornings and nights and less pronounced during the afternoons.

Fig. 4

Changes in apparent hospital position over 10 days due to changes in average refractive index gradient along the path. Observation period is from July 24, 2014 to August 3, 2014.

One interesting aspect of the stitched plot in Fig. 4 is that the hospital appears to be almost at the same observed vertical position every afternoon, suggesting the refractive index gradient values are nearly the same each afternoon throughout the observation period. This can be expected as during the afternoon, the boundary layer is well mixed and the conditions are near adiabatic $(\partial T_{\mathrm{v}}/\partial z=-9.8\mathrm{K}{\mathrm{km}}^{-1}$). A curious dip in the profile beyond the usual floor ($\sim -1\mathrm{m}$) is seen in the early evening of July 27, 2014. In addition to satellite imagery that shows periods of cloudiness on that day, radar and observation data indicate brief periods of rain and a thunderstorm in the afternoon followed by clearing; this led to a strong moisture gradient near the surface. Evaporating rainfall created a virtual temperature near the ground a few Kelvin higher than the dry air temperature; with clearing conditions the moisture near the surface quickly decreased with height as the drier air moved in—this caused the virtual temperature to converge to the dry air temperature aloft and thus decrease with height at a lapse rate slightly less than $-9.8\mathrm{K}{\mathrm{km}}^{-1}$. A “superadiabatic” condition was created, or in other words a layer in the lower atmosphere where the density does not decrease with height as much as it typically does under normal adiabatic conditions. As per Eq. (6), the superadiabatic conditions caused a decrease in curvature for ray paths from the hospital, which caused a drop in the hospital image beyond the level expected for adiabatic conditions. The nighttime excursions in Fig. 4 also show good correlation with meteorological observations. The sharp rise in the profile on the night of July 31 is due to cloud free conditions and weak wind resulting in a very strong thermal inversion. On this night, the temperature increased with height above the ground for the lowest 100 to 200 m (a positive lapse rate). On the contrary, presence of clouds and strong winds did not allow a substantial inversion to develop on the nights of July 26, 27, and 28 which explains the short range of motion seen on those nights. The boundary layer collapses to about 100 m, just above the experimental path in the night. Local disturbances create ripples on top of the stratified boundary layer and this rippling behavior is seen around the nighttime peaks in Fig. 4.

To get an estimate of the refractive index gradient values at different times of the day, observation data for surface temperature and humidity at 2 PM on July 25, 2014, was considered. Assuming the afternoon lapse rate is $-9.8\mathrm{K}{\mathrm{km}}^{-1}$, $n-1=0.000273$ (for 500 nm), and $T_{\mathrm{v}}=295\mathrm{K}$, the vertical refractive index gradient computed using Eq. (6) is $-2.3\times {10}^{-8}{\mathrm{m}}^{-1}$, which is a smaller gradient in the vertical than the standard atmosphere value ($-2.7\times {10}^{-8}{\mathrm{m}}^{-1}$). This is expected, as the lapse rate in the standard atmosphere model is $-6.5\mathrm{K}{\mathrm{km}}^{-1}$, and it implies there is more refractive bending in a standard atmosphere than in a neutrally buoyant atmosphere ($-9.8\mathrm{K}{\mathrm{km}}^{-1}$). Since the gradients were observed to be nearly the same every afternoon, this value of $-2.3\times {10}^{-8}{\mathrm{m}}^{-1}$ was taken as a representative value for the afternoon refractive index gradient during the observation period. Thus, the maximum nighttime values of the average refractive index gradient varied from $-3.5\times {10}^{-8}$ to $-9.5\times {10}^{-8}{\mathrm{m}}^{-1}$.

4.

Conclusions

The U.S. Standard Atmosphere model provides a single value for the atmospheric refractive index gradient. The actual refractive index gradient in any location on any given day varies in correlation with solar heating and other weather conditions. A time-lapse imaging experiment was conducted to estimate the temporal variability of the refractive index gradient. The diurnal variation of the average refractive index gradient was found to be more than twice the standard atmosphere value for this gradient on certain days. The gradient changes were more pronounced during the nights and early mornings. The average refractive index gradient was nearly the same each afternoon throughout the observation period when the atmospheric conditions were nearly adiabatic. Comparisons with observational data showed the gradient variations estimated from the time-lapse imagery correlated well with changing weather conditions too.

The imaging approach can be used as a ground truth measurement of refractive bending for comparison to derived estimation methods, such as those obtained from temperature and moisture gradients. It can thus be used for validation and improvement of fine-scale numerical weather models. Future work will include estimating atmospheric turbulence parameters from the individual and differential random motions of different features in the images.