Retrieving tropospheric refractivity structure using interferometry of aircraft radio transmissions

Detailed measurements of atmospheric humidity in the lower atmosphere are currently difficult and expensive to obtain. For this reason, there is interest in the development of low-cost, high-volume opportunistic technologies to acquire measurements of tropospheric humidity. We demonstrate the use of interferometry to measure the atmospheric refraction of the Automatic Dependent Surveillance-Broadcast (ADS-B) radio transmission routinely broadcast by commercial aircraft. Atmospheric refraction is strongly influenced by changes in humidity, and refractivity observations have proved to be an effective source of humidity information for numerical weather prediction models. A prototype ADS-B interferometer has been developed that can simultaneously perform angle-of-arrival (AoA) interferometry and decode ADS-B signals. Combining the measured AoA of the ADS-B signal with the known position of the aircraft (information contained within the ADS-B signal) allows the bending of the signal due to refraction to be determined. Combining the measured bending of numerous ADS-B signals allows for information concerning the refractivity structure to be extracted. An adjoint model is derived and used to retrieve synthetic one-dimensional refractivity profiles in a variety of atmospheric conditions. The results from an experiment using a prototype ADS-B interferometer are shown, and initial refractivity profiles are retrieved. Sources of uncertainty in the observations and the retrieved refractivity profiles are explored, and future work is suggested.

1 Introduction

The global distribution of water vapour partially governs the transport of heat and the radiation budget of the atmosphere (Trenberth et al., 2009; Held and Soden, 2000). In the lower atmosphere, the spatial and temporal variability in water vapour strongly influences the formation and intensity of extreme weather events through condensation and the subsequent release of latent heat (Trenberth et al., 2005). The high sensitivity of atmospheric processes to water vapour presents the requirement for accurate observations for both climatological and meteorological purposes. Despite its importance in climate and driving extreme weather events, the spatial and temporal variability of tropospheric humidity is difficult to capture adequately using existing observing systems (Rieckh et al., 2018). Radiosondes provide the principal source of in situ tropospheric water vapour observations, capable of retrieving high-resolution vertical profiles of humidity. However, the limited spatial distribution of launch sites and low frequency of radiosonde ascents (typically once daily) substantially restrict the ability to retrieve the full four-dimensional structure of water vapour in the lower atmosphere. Another source of in situ humidity observations is provided by the Aircraft Meteorological Data Relay (AMDAR). The onboard water vapour sensing system (WVSS) instrument can provide detailed profiles of humidity, particularly as the aircraft traverses the moist lower atmosphere (Petersen et al., 2016). However, only a fraction of AMDAR reports currently contain usable humidity data, with the majority restricted to North America (Petersen et al., 2016; Ingleby et al., 2019; Pauley and Ingleby, 2022).

Developments in the remote sensing of atmospheric water vapour have substantially improved the ability of observing systems to capture the variability and distribution of tropospheric humidity. Concerning the remote sensing of humidity, surface-based retrievals using lidars can provide excellent temporal coverage of the variability in humidity throughout the extent of the lower atmosphere directly above the sensor (Guo et al., 2024). However, individual units are expensive and limited in number, restricting the distribution of retrieved vertical profiles of tropospheric humidity for use in numerical weather prediction (NWP) (Gaffard et al., 2021). The maturation of satellite-based remote sensing of humidity has dramatically increased the global coverage of humidity observations. Satellite observations, such as those from passive microwave and infrared-sounding satellites, have proved to be a crucial source of indirect humidity observations. In particular, microwave humidity-sounding has had a substantial positive impact on short-range forecasting skill (Geer et al., 2017). Microwave humidity-sounding satellite observations are an especially crucial source of humidity information over the data-sparse open oceans, where the positive impact of radiance-derived humidity observations in these regions has been demonstrated to propagate downstream with the atmospheric flow over much of Europe and landmasses in the Southern Hemisphere. However, the large weighting functions (2–3 km) used to quantify the vertical distribution of humidity severely limit the resolution of retrieved profiles (Rieckh et al., 2018; Wang et al., 2024) relative to the vertical resolution of NWP models such as the Met Office United Kingdom (UK) Variable-resolution (UKV) model (average vertical spacing in the troposphere is on the order of ∼ 100 m) (Tang et al., 2013).

High-resolution direct/indirect measurements of humidity are challenging to obtain. However, the high spatial and temporal variability of water vapour in the lower atmosphere results in significant changes in the optical properties of the atmosphere. A valuable source of indirect humidity information for use in NWP is provided through observations of refractivity. An example would be the Global Navigation Satellite System (GNSS) radio occultation technique (Kursinski et al., 1997; Kuo et al., 2000; Healy and Thépaut, 2006), which exploits the refractive bending of radio signals between a low Earth orbit (LEO) satellite and a GNSS satellite to retrieve the refractive index (refractivity) structure of the atmosphere. Although GNSS radio occultation can provide accurate high-resolution vertical profiles of refractivity, the retrievals can be significantly impacted by strong gradients in refractivity (Xie et al., 2006). Refractivity observations have also been obtained using radar phase shifts (Fabry et al., 1997; Weckwerth et al., 2005) and zenith tropospheric delay (ZTD) measurements of GNSS satellite radio signals (Rohm, 2012; Zhao et al., 2019). However, the enormous spatial and temporal variability in refractivity makes accurate, high-resolution retrievals challenging.

There is increasing interest in acquiring atmospheric observations using commercial aircraft, which are mandated to broadcast their digitally encoded three-dimensional position via the Automatic Dependent Surveillance-Broadcast (ADS-B) system for air traffic control and safety purposes (ICAO, 2012). Meteorological parameters derived using ADS-B and Mode Selective Enhanced Surveillance (Mode-S EHS) data include temperature (de Haan, 2011; Stone and Kitchen, 2015; Mirza et al., 2019) and wind (de Haan et al., 2013; Mirza et al., 2016) data. Aircraft-derived observations have been exploited by the UK Met Office (Stone and Pearce, 2016; Stone, 2018; Li, 2021) and Royal Netherlands Meteorological Institute (KNMI) (de Haan and Stoffelen, 2012; de Haan et al., 2013) to retrieve high-density wind data. The potential to retrieve atmospheric refractivity data using aircraft-derived observations has recently been explored (Lewis et al., 2023a). A two-element interferometer was used to measure the bending of refracted ADS-B signals, and analogies with GNSS radio occultation were investigated (Lewis et al., 2023b). In this study, we demonstrate the potential to retrieve one-dimensional refractivity profiles using synthetic and real observations obtained using a prototype ADS-B interferometer.

2 Methods

2.1 Interferometry of radio transmissions

Atmospheric refraction of ADS-B radio signals transmitted by aircraft can be detected using vertically orientated, two-element interferometers. The total bending of the radio signal can be determined by combining the interferometrically observed angle-of-arrival (AoA) with the known point of origin of the signals (information encoded within the ADS-B transmission). The observed AoA measurement technique is described in detail in Lewis et al. (2023a). Measurable bending occurs only at grazing incidence (AoA ≲ 2°); therefore, a clear view of the horizon is required to acquire meteorologically useful observations. A measurable phase difference between the signal received at the upper antenna and the signal received at the lower antenna is induced by the path length excess (geometry depicted in Fig. 1). The phase difference, ϕ, can be written in terms of the path excess as

$$\begin{array}{}\text{(1)}& \mathit{\varphi }=\mathrm{2}\mathit{\pi }{\displaystyle \frac{\mathit{B}\cdot \widehat{\mathit{z}}}{\mathit{\lambda }}}=\mathrm{2}\mathit{\pi }{\displaystyle \frac{B\sin \left(\mathit{\beta }\right)}{\mathit{\lambda }}},\end{array}$$

where B is the vertical baseline vector (displacement between the upper and lower antennas), $\widehat{\mathit{z}}$ is the unit vector in the direction of the received ADS-B signal, β is the observed AoA (elevation angle of the received signal) and λ is the signal wavelength (∼ 0.275 m for a 1090 MHz transmission). An Ettus Universal Software Radio Peripheral (USRP) X300 software-defined radio (SDR) was used to determine the phase difference by computing the conjugate product of the signals from each interferometer element.

Figure 1 Geometry of the ADS-B interferometry observing system.

2.2 Atmospheric refraction

The refractive properties of the atmosphere can be described using refractivity, which is the parts-per-million (ppm) excess of the refractive index above the unity vacuum value. The total refractivity at a given point in the atmosphere can be decomposed into the dry refractivity (contribution of the dry atmospheric gases to the total refractivity, N_dry) and the wet refractivity (contribution of water vapour to the total refractivity, N_wet) as

$$\begin{array}{}\text{(2)}& {\displaystyle }{N}_{\mathrm{dry}}=q_{\mathrm{1}}{\displaystyle \frac{P}{T}},\text{(3)}& {\displaystyle }{N}_{\mathrm{wet}}=q_{\mathrm{2}}{\displaystyle \frac{e}{T^{\mathrm{2}}}},\end{array}$$

where q₁ and q₂ are empirically derived constants given by q₁=77.6 K hPa⁻¹ and q₂ = 3.73×10⁵ K² hPa⁻¹, respectively (Smith and Weintraub, 1953). The air pressure, temperature and partial pressure of water vapour are given by P [hPa], T [K] and e [hPa], respectively. Other more accurate formulae exist for the calculation of refractivity (e.g. Aparicio and Laroche, 2011), but the simple form is still commonly used. The contribution of water vapour to the total refractivity is generally around 20 % near the surface (Möller, 2017) but can fluctuate rapidly over short spatial and temporal scales. Refractivity data were obtained using the 12:00 Coordinated Universal Time (UTC) radiosonde launch data from Watnall, Nottinghamshire (UK) (World Meteorological Organisation (WMO) station number 03354) on 22 September 2023. The radiosonde data used in this study are available from the University of Wyoming, United States of America. Figure 2 shows the tropospheric dry refractivity profile, as well as the modulation due to the wet refractivity contribution. Above an altitude of around 7 km, the contribution of the wet refractivity to the total refractivity is negligible.

Figure 2 The total refractivity (solid black line) and dry refractivity (solid green line) versus altitude (in kilometres) retrieved using the 12:00 UTC launch radiosonde data from Watnall on 22 September 2023.

To model the propagation of a radio signal through an atmosphere with a radially directed refractive index gradient, we apply ray-tracing techniques used in geometric optics. The ray path in a spherically symmetric atmosphere can be defined using Snell's law:

$$\begin{array}{}\text{(4)}& n\cdot (a+h)\cdot \cos \left(\mathit{ϵ}\right)=\text{constant},\end{array}$$

where n is the refractive index, h is the height of the ray above the surface, ϵ is the local elevation angle of the ray and a is the radius of the Earth. The ray propagation can be formulated as an initial value problem, dependent on the initial position and direction of the ray. Here, we use the second-order ordinary differential equation (SODE) derived by Zeng et al. (2014) to simulate the propagation. The SODE model is derived using the initial assumption that the height, h, is dependent on the along-beam range, r. Therefore,

$$\begin{array}{}\text{(5)}& n\left[h\left(r\right)\right]\cdot \left(a+h\left(r\right)\right)\cdot \cos \left(\mathit{ϵ}\right)=\text{constant}.\end{array}$$

Differentiating Eq. (5) with respect to r allows the second-order differential equation to be derived (the full derivation is described in detail in Zeng et al., 2014):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}h}{\mathrm{d}{r}^{\mathrm{2}}}}=\left(\mathrm{1}-{\left({\displaystyle \frac{\mathrm{d}h}{\mathrm{d}r}}\right)}^{\mathrm{2}}\right)\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}+{\displaystyle \frac{\mathrm{1}}{a+h}}\right).\end{array}$$

The introduction of the variable u allows the SODE to be separated into a system of first-order differential equations:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}h}{\mathrm{d}r}}=u,\end{array}$$

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}u}{\mathrm{d}r}}=\left(\mathrm{1}-u^{\mathrm{2}}\right)\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}+{\displaystyle \frac{\mathrm{1}}{a+h}}\right),\end{array}$$

which can be integrated over the domain [0,r_f] (where r_f defines the endpoint of the ray) to recover the ray path. We choose to integrate Eqs. (7) and (8) using a third-order Runge–Kutta scheme. The distance travelled by the ray at each step, i, along the surface of the Earth, s, is given by

$$\begin{array}{}\text{(9)}& s_{i+\mathrm{1}}=s_i+{\sin }^{-\mathrm{1}}\left({\displaystyle \frac{\cos \left({\mathit{ϵ}}_i\right)\mathrm{\Delta }r}{a+h_{i+\mathrm{1}}}}\right)a,\end{array}$$

where

$$\begin{array}{}\text{(10)}& {\mathit{ϵ}}_i={\sin }^{-\mathrm{1}}\left(u_i\right)\end{array}$$

is the local elevation angle of the ray at each step. The initial height and direction of the ray are required to model the refraction of a radio signal through the atmosphere. The altitude of the aircraft can be determined; however, the initial emission angle of the radio signal is unknown. The refraction is therefore modelled in a time-reversed frame from the position of the receiver to the reporting aircraft. A prior (first-guess) refractivity profile is used. The squared difference between the altitude of the model ray after termination and the reported height of the aircraft is used to compute a penalty function. The best estimate of the true refractivity profile is determined by minimising the penalty function. We explore this optimisation problem in more detail in Sect. 2.3.

2.3 Optimisation of the refractive index profile

The optimisation problem can be formulated as

$$\begin{array}{}\text{(11a)}& {\displaystyle }\underset{n}{min}J=\sum _{i=\mathrm{0}}^{N_{\mathrm{B}}}{j}_i=\sum _{i=\mathrm{0}}^{N_{\mathrm{B}}}\left(\right(h_{\mathrm{f}}{)}_i-(h_{\mathrm{t}}{)}_i{)}^{\mathrm{2}},\text{(11b)}& {\displaystyle }\mathrm{\text{subject to}}\dot{h}(r;n,h_{\mathrm{0}},u_{\mathrm{0}})-u=\mathrm{0},\forall r\in [\mathrm{0},r_{\mathrm{f}}],\text{(11c)}& {\displaystyle }\begin{array}{l}\dot{u}(r;n,h_{\mathrm{0}},u_{\mathrm{0}})-(\mathrm{1}-u^{\mathrm{2}})\left(\frac{\mathrm{1}}{n}\frac{\mathrm{d}n}{\mathrm{d}h}+\frac{\mathrm{1}}{a+h}\right)=\mathrm{0},\\ \forall r\in [\mathrm{0},r_{\mathrm{f}}],\end{array}\text{(11d)}& {\displaystyle }h(\mathrm{0};n,h_{\mathrm{0}},u_{\mathrm{0}})=h_{\mathrm{0}},\text{(11e)}& {\displaystyle }u(\mathrm{0};n,h_{\mathrm{0}},u_{\mathrm{0}})=u_{\mathrm{0}},\end{array}$$

where h and u describe the height (position) and (sine) direction of the ray, respectively, and $\dot{h}$ and $\dot{u}$ are their respective derivatives with respect to r. The initial position and direction of the ray are given by h₀ and β, respectively. J is the total penalty function describing the sum of the individual penalty functions, j_i, corresponding to the squared difference between the end height of the ith ray, (h_f)_i, and the ith target height, (h_t)_i ($i=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{B}}$, where N_B is the number of broadcasts used in the optimisation problem). In the case of ADS-B interferometry, the target position will be the reported position of the aircraft.

Figure 3 shows the geometry of the system to be optimised for a single ray. The refractive index profile, n, was optimised such that a ray traced through the profile with an initial position and direction given by h₀ and β (where u₀ = sin (β)), respectively, will terminate at the target (aircraft) position, h_t. The gradient of the penalty function, J, with respect to the refractive index profile, n, was required to minimise J via gradient descent. The gradient of the individual penalty function for a single ray can be approximated by computing finite differences over the elements of the refractive index profile n∈ℝ^M (where M is the number of refractive index values along the ray trajectory) as

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{d}j\left[h_{\mathrm{f}}\left(\mathit{n}\right)\right]}{\mathrm{d}{n}_k}}\approx {\displaystyle \frac{j\left[h_{\mathrm{f}}(\mathit{n}+\mathit{\delta }n{\mathit{e}}_k)\right]-j\left[h_{\mathrm{f}}(\mathit{n}-\mathit{\delta }n{\mathit{e}}_k)\right]}{\mathrm{2}\mathit{\delta }{n}_k}},\end{array}$$

where n_k is the kth element of the refractive index profile. The kth element of the vector e_k (e_k∈ℝ^M) is unity and zero otherwise, and δn is a small increment to the kth element of n. Each gradient evaluation with respect to a single element of the refractive index profile requires integrating Eqs. (11b) and (11c) over the entire ray trajectory to determine the change in the end position of the ray, h_f. The penalty gradient calculation using finite differences for a single ray is expensive and quickly becomes impractical when evaluating the total penalty function gradient for N_B rays when N_B is large.

In a previous investigation by Lewis et al. (2023a), the total penalty function given by Eq. (11a) was minimised using simulated annealing, a stochastic global search optimisation algorithm (Kirkpatrick et al., 1983). Although simulated annealing is an effective method for optimising high-dimensional problems with potentially many local minima, the algorithm is very computationally expensive and requires thousands of iterations to converge on a good solution. For this reason, the optimisation procedure using the adjoint state method (Zou et al., 1997; Gill et al., 1981) was used to reduce the computational cost of retrieving the refractivity profile. Applying the adjoint state method to our original optimisation problem, the constrained optimisation problem described by Eq. (11) can be converted into an unconstrained optimisation problem using Lagrange multipliers:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\mathcal{L}& =j-\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\psi }\left(\dot{h}-u\right)\mathrm{d}r\\ & -\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\mu }\left(\dot{u}-(\mathrm{1}-u^{\mathrm{2}})\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}+{\displaystyle \frac{\mathrm{1}}{a+h}}\right)\right)\mathrm{d}r,\end{array}\end{array}$$

where ψ and μ are the adjoint variables associated with the state variables h and u, respectively. The integration of equations Eqs. (11b) and (11c) over the domain [0,r_f] (where r_f parameterises the endpoint of the ray) defines the trajectory of the ray. The minimisation of Eq. (13) corresponds to finding the stationary points of the Lagrangian with respect to the trajectory defined by h and u. Therefore, any variations in the Lagrangian with respect to the trajectory are equal to 0 (that is, δℒ|_h = δℒ|_u = 0). The adjoint variables enforce this constraint. Solving for the stationary points of the Lagrangian with respect to h and u, we find

$$\begin{array}{}\text{(14)}& \mathit{\psi }{|}_{r=r_{\mathrm{f}}}={\displaystyle \frac{\partial j}{\partial h}},\end{array}$$

$$\begin{array}{}\text{(15)}& \mathit{\mu }{|}_{r=r_{\mathrm{f}}}={\displaystyle \frac{\partial j}{\partial u}},\end{array}$$

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}r}}=\mathit{\mu }(\mathrm{1}-u^{\mathrm{2}})\left({\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}\right)}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}n}{\mathrm{d}{h}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{(a+h{)}^{\mathrm{2}}}}\right),\end{array}$$

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}r}}=\mathrm{2}\mathit{\mu }u\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}+{\displaystyle \frac{\mathrm{1}}{a+h}}\right)-\mathit{\psi }.\end{array}$$

The full derivation is described in Appendix A. The initial conditions of the adjoint variables ψ and μ are defined at the endpoint of the ray trajectory, $h{|}_{r=r_{\mathrm{f}}}$. The adjoint equations are therefore evolved in reverse from the endpoint to the initial position of the ray. In the simulated retrievals, an atmosphere with spherical symmetry was assumed, and values were computed using a linear interpolation scheme; hence, ${\mathrm{d}}^{\mathrm{2}}n/\mathrm{d}{h}^{\mathrm{2}}$ = 0. Because the refractive index generally decreases exponentially with height, the interpolation was performed using the natural logarithm of the refractive index values. Expressing the height-dependent natural logarithm of the refractive index as $\stackrel{\mathrm{̃}}{n}\left(h\right)$ = ln (n(h)), Eqs. (16) and (17) can then be written as

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}r}}=\mathit{\mu }(\mathrm{1}-u^{\mathrm{2}})\left({\left({\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{̃}}{n}}{\mathrm{d}h}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{(a+h{)}^{\mathrm{2}}}}\right),\end{array}$$

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}r}}=\mathrm{2}\mathit{\mu }u\left({\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{̃}}{n}}{\mathrm{d}h}}+{\displaystyle \frac{\mathrm{1}}{a+h}}\right)-\mathit{\psi }.\end{array}$$

Figure 3 Illustration of the single-ray optimisation problem. A ray with initial position, h₀, and elevation (exaggerated in the diagram for illustrative purposes), β, is traced through an initial-guess refractive index profile until it reaches the aircraft distance, s_t. The penalty function for a single ray is calculated as the squared difference between the ray endpoint, h_f, and the target (aircraft) position, h_t. The intermediate height, distance and local elevation angle of the ray along the path are given by h_i, s_i and ϵ_i, respectively.

The integration of the adjoint variables requires knowing the spatial gradient of the refractive index at the query points h along the ray trajectory. The refractive index profile was optimised using the two-stage optimisation procedure described by Teh et al. (2022): forward tracing and backward tracing. The forward tracing stage involves integrating Eqs. (11b) and (11c) using initial conditions h₀ and u₀ and evaluating the endpoint of the ray trajectory $h{|}_{r=r_{\mathrm{f}}}$. The backward tracing stage then involves initialising the adjoint variables given by Eqs. (14) and (15) at the termination position and integrating Eqs. (18) and (19) in reverse along the same ray trajectory as in the forward tracing stage. The variation of the Lagrangian with respect to the refractive index is then

$$\begin{array}{}\text{(20)}& \mathit{\delta }\mathcal{L}{|}_{\stackrel{\mathrm{̃}}{n}}=\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\mu }(\mathrm{1}-u^{\mathrm{2}})\left({\displaystyle \frac{\mathrm{d}\left(\mathit{\delta }\stackrel{\mathrm{̃}}{n}\right)}{\mathrm{d}h}}\right)\mathrm{d}r,\end{array}$$

which does not require an expensive Jacobian matrix to evaluate. The quantity $\mathit{\delta }\stackrel{\mathrm{̃}}{n}$ is defined in terms of a linear interpolation scheme as

$$\begin{array}{}\text{(21)}& \stackrel{\mathrm{̃}}{n}(h,\mathit{n})=\mathit{w}\left(h\right)\cdot \ln \left(\mathit{n}\right)=\left[\begin{array}{cc}{w}_{\mathrm{0}}\left(h\right)& w_{\mathrm{1}}\left(h\right)\end{array}\right]\cdot \left[\begin{array}{c}\ln \left(n_{\mathrm{0}}\right)\\ \ln \left(n_{\mathrm{1}}\right)\end{array}\right],\end{array}$$

where n₀ and n₁ are the refractive index values on the underlying grid directly below and above the sample point h, respectively. The weights w₀(h) and w₁(h) are inversely related to the vertical distance from the refractive index values n₀ and n₁, respectively. Therefore, at each step in the integration of the adjoint variables, the gradient of the Lagrangian with respect to the refractive index values on the underlying grid is calculated as

$$\begin{array}{}\text{(22a)}& {\displaystyle \frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\left(\ln \left(n_{\mathrm{0}}\right)\right)}}=\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\mu }\left(\mathrm{1}-u^{\mathrm{2}}\right){\displaystyle \frac{\mathrm{d}{w}_{\mathrm{0}}}{\mathrm{d}h}}\mathrm{d}r,\text{(22b)}& {\displaystyle \frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\left(\ln \left(n_{\mathrm{1}}\right)\right)}}=\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\mu }\left(\mathrm{1}-u^{\mathrm{2}}\right){\displaystyle \frac{\mathrm{d}{w}_{\mathrm{1}}}{\mathrm{d}h}}\mathrm{d}r,\end{array}$$

where the relation $\mathit{\delta }\stackrel{\mathrm{̃}}{n}$ = w⋅d(ln(n)) is employed. The geometry of the optimisation scheme at each step in the integration is shown in Fig. 4.

Figure 4 The geometric ray path intersecting the refractive index values ln (n₀) and ln (n₁) on the underlying grid.

The refractive index field was adjusted iteratively using a gradient descent scheme using a variation of the Adam optimiser (Kingma and Ba, 2014). The minimisation algorithm was terminated after a specified number of iteration steps. The gradient of the penalty function calculated using the adjoint state method was verified using finite differences over the elements of the refractive index profile, as shown in Fig. 5. The penalty function is highly sensitive to the surface refractive index value (as shown in Fig. 5a, c), highlighting the importance of accurately knowing the surface conditions in order to constrain the optimised solution. Figure 5b, d show the magnitude of the relative difference between the gradients calculated using finite differences and the adjoint state method for 1 and 100 m ray integration step sizes, respectively. The value of the gradients computed using finite differences and the adjoint state method converge as the ray step size decreases; however, the computational cost of integrating rays with a very small step size increases significantly. For this study, the 100 m integration step size was used, as the gradient calculation was sufficiently accurate without being prohibitively expensive to compute for thousands of rays.

Figure 5 (a, c) The gradient of the penalty function with respect to the refractive index at each altitude point calculated using finite differences (green) and the adjoint state method (black) for a (a) 1 m and (c) 100 m ray integration step size. The penalty function is highly sensitive to the refractive index values close to the surface. (b, d) The magnitude of the relative difference between the gradients calculated using finite differences (FD) and the adjoint state method (adj) for a (b) 1 m and (d) 100 m ray integration step size. The ray had an observed AoA of 0.2°. The altitude and distance of the source aircraft were 8.8 and 350 km, respectively.

In summary, there are two principal observations required to retrieve the atmospheric refractivity structure: the observed AoA of the radio broadcast (required to calculate the initial ray direction, u₀) and the position of the aircraft (required to determine s_t and h_t). In the examples above, we initialise the adjoint variable μ using the partial derivative of the penalty function with respect to the ray direction u. However, in practice, we have no knowledge of the emission direction of the ADS-B signal. Therefore, we define the penalty function only in terms of the respective positions of the ray endpoint and the aircraft position. The adjoint term μ is therefore initialised at the termination position as $\mathit{\mu }{|}_{r=r_{\mathrm{f}}}=\mathrm{0}$.

2.4 Experimental arrangement

The observations were obtained using a prototype ADS-B interferometer installed on the Clee Hill weather radar tower in Shropshire, UK. The experimental setup is shown in Fig. 6. The vertical baseline of the interferometer was 13.86 m, and each of the two elements consisted of a stack of four Yagi–Uda antennas separated by 0.66 m. The stack reduced the beamwidth of each element, helping to mitigate the effects of signals reflected from nearby topography (multipath). The antenna stacks were connected to an Ettus USRP X300 SDR.

Figure 6 The interferometer mounted to the side of the Clee Hill weather radar tower.

Approximately 263 500 refracted ADS-B transmissions were received between an AoA of 0.0 to 2.0° over an ∼ 80 min period between 08:36 and 09:58 UTC on 22 September 2023. The synthetic and real retrievals were obtained using the central 10° azimuthal sector of received transmissions. The distribution of all the received transmissions and those used for retrievals are shown in Fig. 7a and b, respectively. The distribution of reported aircraft altitudes within the chosen sector is shown in Fig. 7c. The northeast sector was chosen due to a suitable distribution in reported altitudes for the lower atmosphere to be adequately sampled. The nearby Watnall radiosonde station also allowed for a reference refractivity ground truth to be determined to compare with retrieved refractivity profiles.

Figure 7 (a) The spatial distribution of all the received aircraft transmissions colour-coded by the reported altitude (km). (b) The spatial distribution of the analysed aircraft transmissions colour-coded by the reported altitude (km). The positions of the interferometer and the Watnall radiosonde station are indicated by the black star and black circle, respectively. (c) Histogram of the distribution of the analysed broadcasts by the reported altitude (km). The analysed aircraft transmissions were sampled from the central 10° azimuthal sector of received transmissions on 22 September 2023.

The adjoint inversion algorithm was tested on synthetic data generated with ADS-B data recorded using the prototype interferometer at Clee Hill, Shropshire. The observations used in the analysis were restricted to the central 10° azimuthal sector of the entire observation dataset. A subset of 5000 ADS-B transmissions were chosen randomly from within the central azimuthal sector to use in the synthetic refractivity retrieval experiments. Radiosonde data from the nearby Watnall station in Nottinghamshire (position shown in Fig. 7a) were used to generate a refractivity profile to simulate the refraction of the subset of ADS-B transmissions. Three radiosonde soundings were used to generate the synthetic observations – 12:00 UTC on 22 September 2023 (see Fig. 2), 12:00 UTC on 18 July 2022 (during an intense heatwave) and 12:00 UTC on 15 December 2022 (during a cold spell). The synthetic observations (heights of synthetic aircraft) were generated by tracing the 5000 rays out from the receiver to the reported distances, s_t (see Fig. 3), of each of the ADS-B transmissions using the refractivity profiles determined using the radiosonde data. The interferometrically measured observed AoAs were used as the initial directions of the rays, and the height of the ray endpoint was used as the “true” height of the synthetic ray origin (i.e. equivalent to the height of a synthetic aircraft). The method is equivalent to that shown in Fig. 3, except the heights of the synthetic aircraft are given by h_f. The heights of the synthetic aircraft are treated as the target heights h_t in the synthetic retrievals.

The observed AoAs used in the synthetic refractivity retrievals were injected with varying degrees of noise to model the impact of experimental uncertainties on the retrieved refractivity profiles. A prior (first-guess) refractivity profile, given by

$$\begin{array}{}\text{(23)}& N\left(h\right)=N_{\mathrm{0}}\exp \left(-{\displaystyle \frac{h}{H}}\right),\end{array}$$

was used to initialise the penalty function, where N₀ is the refractivity at the position of the receiver (this was fixed to the radiosonde value during the optimisation procedure) and H is the scale height in kilometres. Rays were traced backwards from the receiver position to the reported distance of the aircraft using the known observed AoAs, as shown in Fig. 3. The total penalty function was then evaluated by summing together the penalty function for each ray (i.e. the squared difference between the traced ray end height and the height of the synthetic aircraft). Because multipath contamination of the observed AoA measurement is thought to constitute the principal source of uncertainty in the retrieved refractivity profile, three observed AoA measurement noise test cases were investigated. The retrieval algorithm was run with a zero measurement noise case and two normally distributed noise test cases (with standard deviations of 0.01 and 0.05°, respectively).

To account for the ellipsoid shape of the Earth in the one-dimensional ray tracing model, the mean local radius of the ellipsoid at the location of the receiver was used as the effective radius of the Earth. The effective radius is approximately the radius of curvature of the ellipsoid at the location of the receiver, which ensures that contours of constant refractivity are approximately parallel to the curved surface of the ellipsoid. The height above the ellipsoid reported in the ADS-B transmission was left unchanged. The coordinate transform is described in Appendix B.

3 Results and analysis

3.1 Synthetic refractivity retrievals

The retrieved refractivity profile consisted of 30 levels on a logarithmically spaced grid between an altitude of 0.575 km (height of receiver above the ellipsoid Earth model, estimated using a GNSS receiver on site) and an altitude of ∼ 13 km. A ray step size of Δr = 0.1 km was used during the integration. A learning rate of 10⁻⁷ was used during the optimisation procedure. In order to better constrain the retrieved solution (particularly in the presence of noise), the minimum refractivity after every iteration during the optimisation procedure was fixed to be greater than or equal to the dry refractivity (i.e. assuming pressure and temperature are known and the relative humidity can never drop below 0 %). The root mean squared error (RMSE) values in the synthetic retrieved refractivity profiles after convergence for the test cases using measurement noises of 0.00, 0.01 and 0.05° are shown in Table 1. The RMSE was computed as

$$\begin{array}{}\text{(24)}& \text{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{\mathrm{30}}{\left(N_{\mathrm{retrieve},i}-N_{\mathrm{radiosonde},i}\right)}^{\mathrm{2}}}{\mathrm{30}}}},\end{array}$$

where N_retrieve is the retrieved refractivity and N_radiosonde is the refractivity determined from radiosonde data.

Table 1 The RMSE values for the initial refractivity profile and the synthetic retrieved refractivity profile (relative to the reference radiosonde profile). The three angular values refer to the standard deviation of the added noise to the observed AoA.

The RMSE in the retrieved profile increases as the standard deviation in the observed AoA noise increases. However, all retrieved profiles show improvement relative to the prior refractivity profile. A synthetic retrieval with no observed AoA measurement noise was also performed assuming the height of the receiver above the ellipsoid had an error of ±5 and ±10 m. The synthetic retrieval assuming a receiver height error still showed significant improvement in the lower relative to the initial refractivity profile. However, there was little improvement above an altitude of ∼ 10 km. The weakly constrained refractivity retrievals at high altitudes are explored in Sect. 4.

Figure 8a shows the synthetic retrieved vertical refractivity profile for the different observed AoA measurement noise cases using the 22 September radiosonde data. The 0.00° RMSE profile was obtained from the 5000 interferometrically measured, unperturbed observed AoAs. The 0.01 and the 0.05° RMSE profiles were obtained by perturbing the 5000 interferometrically measured observed AoAs according to the corresponding normal distributions. The initial refractivity profile was defined using Eq. (23) with a scale height H = 8 km. The percentage difference between the retrieved and the radiosonde-derived refractivity profile is shown in Fig. 8b. The 0.05° RMSE profile shows a lower refractivity error standard deviation at higher altitudes than the 0.01° RMSE profile. The increased standard deviation in the AoA measurement error results in a greater number of AoAs being perturbed below 0° for the 0.05° RMSE case. Because the retrievals can currently only be performed using AoAs ≥ 0°, these are then rejected, resulting in a bias towards perturbing the AoAs towards higher elevation angles. This results in the gradient in the retrieved refractivity profile becoming steeper, which, as a consequence of the clamping of minimum refractivity, results in the 0.05° RMSE profiles clustering together with a lower variance. However, the larger RMSE in the observed AoA results in a larger bias in the retrieved refractivity profile.

Figure 8 (a) The synthetic retrieved refractivity profile and the reference 22 September 2023 Watnall radiosonde profile for different synthetic error distributions injected into the observed AoA measurement. The error bars indicate the range of retrieved profiles using 20 different observed AoA noise realisations. (b) The percentage difference between the synthetic retrieved refractivity profile and the reference radiosonde profile for different synthetic error distributions injected into the observed AoA measurement. The shaded regions indicate the limits of uncertainty in the synthetic retrieved refractivity profile assuming no AoA measurement noise and a ± 5 m (red) and ± 10 m (blue) receiver altitude error. The fractional error standard deviation in the retrieved refractivity profiles using 20 different observed AoA noise realisations is given by the blue and red dashed lines for the noise distributions with standard deviations of 0.01 and 0.05°, respectively.

The synthetic refractivity retrievals were then repeated for two other radiosonde profiles – the 12:00 UTC Watnall soundings for 18 July 2022 and 15 December 2022. The synthetic retrievals were performed using synthetic observations of aircraft heights generated from the two respective refractivity profiles in a process identical to that using the 22 September data. The synthetic retrievals using injected measurement noise with a standard deviation of 0.01 and 0.05° both showed improvements relative to the initial guess refractivity profile in the lower atmosphere for the 18 July and 15 December profiles, as shown in Table 1. The synthetic refractivity retrievals demonstrate the ability to obtain useful atmospheric data using measurements of refracted ADS-B transmissions in the presence of random noise.

3.2 Initial retrievals using observations

The retrieval algorithm was then tested on the observational data collected using the prototype ADS-B interferometer. The retrieved refractivity profiles were constrained by fixing the minimum refractivity to the dry refractivity. It is hoped that the technique will be able to resolve variations in humidity structure occurring on timescales on the order of an hour. Accordingly, the observed AoA measurements from 1 h of the experiment were divided into four 15 min periods between 08:36 and 09:36 UTC and used to retrieve the refractivity profiles. The observed AoA ranged between 0.0 and 2.0°, and the azimuthal angle range was chosen as the central 10° of the sector shown in Fig. 7. The four observational periods contained 9700 (08:36–08:51 UTC), 8315 (08:51–09:06 UTC), 8185 (09:06–09:21 UTC) and 7589 (09:21–09:36 UTC) received ADS-B transmissions within the sector. Figure 9 shows the

• AoA_LoS: the line-of-sight (LoS) AoA determined from the reported position of the aircraft;

• AoA_{init_model}: the synthetic LoS AoA determined from the end position of the ray traced through the initial refractivity profile;

• AoA_{ret_model}: the synthetic LoS AoA determined from the end position of the ray traced through the retrieved refractivity profile;

as a function of azimuthal angle for the four 15 min periods of observational data. The LoS AoA refers to the angle above the local horizon of the straight line joining the receiver and the reported aircraft position/end position of a traced ray (the LoS AoA computation is explained in more detail in Lewis et al., 2023a).

Figure 9 The LoS AoA as a function of azimuthal angle (where 0° corresponds to the centre of the observational sector). The black points (AoA_LoS) indicate the LoS AoA derived from the decoded ADS-B transmissions. The blue points indicate the modelled LoS AoA (AoA_{init_model}) determined from tracing rays through the initial (first-guess) refractivity profile. The green points indicate the modelled LoS AoA (AoA_{ret_model}) determined from tracing rays through the retrieved (optimised) refractivity profile. The time period was between 08:36 and 09:36 UTC.

The penalty function is minimised when the difference between the AoA_LoS and the AoA_{ret_model} is minimised. The AoA_{init_model} and AoA_{ret_model} simulated using the model refractivity profiles (as shown by the blue and green points in Fig. 9, respectively) provide a direct comparison with observable quantities (the AoA_LoS determined using the received ADS-B transmission). The AoA_{init_model} and AoA_{ret_model} show a ripple-like structure due to significant multipath. To model the LoS AoA using a given refractivity profile, the ray is traced using the interferometrically measured observed AoA (as shown in Fig. 3). Because the observed AoA is noisy, the modelled LoS AoA will share a similar noisy structure. The impact of multipath on the retrieved refractivity profiles is explored below. However, despite the presence of noise in the observed AoA measurement, the AoA_{ret_model} generally converges towards the AoA_LoS. The corresponding refractivity profiles for each 15 min sample of observations between 08:36 and 09:36 UTC are shown in Fig. 10. The reference refractivity profile was determined using the 12:00 UTC Watnall radiosonde launch on 22 September 2023. The refractivity value at the altitude of the receiver was fixed to the radiosonde value at an altitude of 0.575 km. The refractivity profiles retrieved using observational data show agreement with the reference profile radiosonde profile in the lower atmosphere up to an altitude of ∼ 5 km. At higher altitudes, the refractivity profiles show significant variability compared to lower altitudes, likely due to a lack of constraints here (a change in refractivity at higher altitudes has less impact than a change in refractivity at lower altitudes). The potential systematic uncertainties in the observations suggested by these results are explored in Sect. 4.

Figure 10 The 12:00 UTC 22 September 2023 Watnall radiosonde refractivity profile (dotted) and the retrieved (solid colour) and initial (dashed) refractivity profiles. The time period was between 08:36 and 09:36 UTC. The four observational periods contained 9700 (08:36–08:51 UTC), 8315 (08:51–09:06 UTC), 8185 (09:06–09:21 UTC) and 7589 (09:21–09:36 UTC) received ADS-B transmissions within the sector.

The difference between the AoA_LoS and AoA_{init_model} and the difference between the AoA_LoS and AoA_{ret_model} as a function of surface distance for the four observational periods are shown in Fig. 11. The distribution for the difference between the AoA_LoS and AoA_{init_model} and the difference between the AoA_LoS and AoA_{ret_model} is shown in Fig. 12. The mean of the LoS AoA difference decreased from ∼ 0.034 to ∼ 0.009° and the standard deviation increased from ∼ 0.023 to ∼ 0.024° before and after the retrieval, respectively. The standard deviation in the LoS AoA difference exceeds the target accuracy of 0.01°, meaning that a complex refractivity structure is unlikely to be resolved. The AoA_{init_model} is larger than the AoA_LoS, indicating that the initial atmosphere is less refractive than the true atmosphere. The retrieved atmosphere results in apparent weaker refraction of ADS-B signals from nearby aircraft and apparent stronger refraction for signals originating from more distant aircraft. The apparent variability in the observed refraction of ADS-B signals originating from aircraft at varying distances is likely to be partially a consequence of observational noise. The noise is visible in Fig. 9 as a ripple-like structure in the flight paths. This noise is correlated with the foreground environment and is not truly random, which likely results in a bias in the retrieved refractivity profiles. Potential sources of uncertainty, including multipath interference and instrumental errors, are described in Sect. 4.

Figure 11 The difference between the AoA_LoS and AoA_{init_model} (blue line) and the difference between the AoA_LoS and AoA_{ret_model} (green line) as a function of distance. The time period was between 08:36 and 09:36 UTC.

Figure 12 The distribution for the modelled LoS AoA difference for all observations using the initial (blue) and retrieved (green) refractivity profiles. The mean difference in the LoS AoA decreased from ∼ 0.034 to ∼ 0.009° for the initial and retrieved refractivity profiles, respectively. The standard deviation in the LoS AoA difference increased from ∼ 0.023 to ∼ 0.024° for the initial and retrieved refractivity profiles, respectively.

4 Discussion

Initial tests with synthetic data and the case study using actual measurements have revealed possible weaknesses in the retrieval method. Sources of experimental uncertainty that require further attention are discussed in the subsections below.

4.1 Assumption of a horizontally homogeneous refractivity field

The SODE ray tracing model was derived with the assumption of spherical symmetry (the refractivity varies only with radial distance and is horizontally homogeneous). Therefore, only one-dimensional refractivity profiles can be derived using the adjoint model presented here. For rays with negative elevation angles, this assumption will likely break down rapidly as the ray propagates through the highly variable planetary boundary layer. In this study, only ADS-B transmissions with an observed AoA ≥ 0° were used. As the rays generally increase in altitude with distance due to the curvature of the Earth, the assumption of spherical symmetry is therefore assumed to be reasonable. However, significant horizontal variability in refractivity can occur along frontal systems, which can invalidate the spherical symmetry assumption. An occluded front moved southeastwards during the observation period, which may have introduced some stronger horizontal variability than usual. The weather during the observation period was partially cloudy, with patchy low cloud present around Clee Hill. Local variations in the refractivity near the receiver could contribute to uncertainty in the refractivity retrievals. However, the significant multipath contamination of the observed AoA measurement was expected to be the dominant source of uncertainty in the retrieved refractivity profiles.

To assimilate ADS-B refraction measurements influenced by significant horizontal refractivity gradients, a two-dimensional ray tracing model and its adjoint will be required. It is uncertain how well constrained a two-dimensional refractivity retrieval will be using a single ADS-B interferometer. A network of ADS-B interferometers across the UK could potentially allow for two- or three-dimensional refractivity profiles to be obtained.

4.2 Sensitivity of the retrievals to surface refractivity

The retrieved refractivity profile is extremely sensitive to the surface refractivity (as shown in Fig. 5). Uncertainty in the measured refractivity here will result in erroneous retrievals; therefore, an accurate in situ measurement of refractivity will likely be required to constrain the retrieved solution. Accurate on-site refractivity measurements could be obtained using a permanent weather station. The bending of radio signals is dominated by the refractivity gradient rather than the absolute refractivity (as the refractive index of air is approximately unity even near the surface). Any uncertainty in the surface refractivity introduces a constant offset throughout the entire retrieved profile; however, the retrieved refractivity gradient remains unchanged. The surface refractivity uncertainty could be addressed by instead assimilating refractivity gradients into NWP models, which would reduce the need for accurate surface refractivity measurements.

4.3 Multipath

Estimating the observed AoA to an uncertainty of 0.01° requires accurate measurements of the phase difference across the interferometer elements. A significant source of uncertainty in this measurement is due to the presence of multipath. At grazing elevation (observed AoA ≲ 2.0°), the direct and reflected ADS-B signals are received in rapid succession, resulting in mutual interference. The multipath results in ripple-like structures in the observed flight paths of the aircraft, as shown in Fig. 9. The uncertainties in the observed AoA measurements are highly correlated and dependent on the local terrain, which makes estimating the observation error covariance matrix very difficult. Efforts to minimise the impacts of multipath have included beam forming by stacking multiple antennas on each element of the interferometer (see Fig. 6) and increasing the sampling rate. The beam forming technique aims to deepen and control the position of the antenna nulls to minimise the sensitivity to foreground reflections. Increasing the sampling rate allows for phase measurements to be extracted from the direct signal before the arrival of the reflected signal. The sampling rate used for the experiment was 10 MHz, which was able to strongly suppress multipath originating from foreground reflections with a delay longer than ∼ 100 ns (excess path length ∼ 30 m). Multipath due to foreground reflections with a path excess shorter than this can only be partially suppressed using the current ADS-B interferometer. Future experiments will aim to further suppress the multipath by increasing the sampling rate and modifying the beam pattern of the stacks. Both techniques have resulted in significant improvements in the accuracy of the observed AoA measurement. The multipath due to foreground reflections could also be mitigated using a long time series of observations. The impact of reflections from a static foreground would be a systematic perturbation to the observed phase difference. With enough observations from aircraft traversing at a range of elevations and distances, the constant multipath interference could be well characterised after accounting for atmospheric variability. The non-constant variability, such as atmospheric refractivity changes, could be detrended out with sufficient observations. A similar technique is used in ground-based GNSS radio occultation, where Padullés et al. (2016) describes how time series of the observations can be used to remove a significant contribution of the multipath effect on the uncertainty in the measurements of the atmospheric polarimetric effects on observables. Future experiments using the prototype ADS-B interferometer will aim to quantify the nature of the observational noise arising due to multipath contamination.

4.4 Instrumental uncertainties

Instrumental errors may include the assumption that a measured phase difference of 0 corresponds to an observed AoA of 0.00°. A slight offset between the upper and lower elements of the interferometer could introduce uncertainty in this assumption. The reported position of the aircraft is also limited by the reporting resolution (on the order of ∼ 10 m; Stone and Kitchen, 2015); however, the impact of positional uncertainty on the observed AoA was expected to be minor due to the large distances to the aircraft (Lewis et al., 2023a). Uncertainty in the altitude of the receiver will introduce a systematic error in the retrieved refractivity profiles. Figure 8 shows the impact of uncertainty in the height of the receiver on the retrieved refractivity profiles. The height uncertainty results in significant variability in refractivity at high altitudes. The weakly constrained refractivity retrievals above an altitude of ∼ 5 km for some of the retrieved profiles in Fig. 10 suggest significant observational noise. This suggests a potential systematic error in the observational data. However, because the high altitude variability was only present in some of the retrieved profiles, it may suggest that the source of the uncertainty may not be entirely systematic. The gradient of the penalty function with respect to refractivity at higher altitudes is significantly lower than at lower altitudes, as shown in Fig. 5. Introducing additional constraints to the penalty function to penalise large departures from the reference refractivity profile may reduce the presence of unphysical refractivity structures. In practice, a one-dimensional variational (1D-Var) data assimilation scheme would prevent unphysical retrievals. However, a good understanding of the observational uncertainties would be required.

The mean altitude of the midpoint between the elements of the interferometer above the ellipsoid Earth model was estimated to be 575 m using a GNSS receiver located on-site. The standard deviation of the measured altitude was estimated to be 2 m. Multipath can result in correlated uncertainties in the observed AoA measurements, potentially resulting in biased retrievals.

5 Conclusions

A new method has been developed to retrieve refractivity profiles from the interferometrically measured AoA of radio signals from aircraft. The method incorporates the SODE ray tracing model (Zeng et al., 2014) and was tested on both synthesised data (to explore sensitivity to measurement noise) and real observational data measured with an experimental interferometer. Three AoA measurement noise cases were simulated, consisting of no noise and two normally distributed noise cases (with standard deviations of 0.01 and 0.05°, respectively) injected into the observed AoA measurements. The impact of each of the injected observed AoA measurement uncertainty cases on the retrieved vertical refractivity profiles was investigated. The convergence of the penalty function was impacted negatively, and the RMSE in the retrieved refractivity profile increased as the standard deviation in the injected (normally distributed) noise increased. The synthetic retrievals using an AoA noise standard deviation of 0.01° were of comparable accuracy to vertical refractivity profile retrievals using GNSS radio occultation sounding (Anthes et al., 2022), airborne GNSS radio occultation sounding (Murphy et al., 2015) and GNSS troposphere tomography (Wu et al., 2014; Trzcina and Rohm, 2019). When the algorithm was tested on real measurements in four separate 15 min periods, the retrieved profiles showed significant variability at high altitudes but good convergence towards the reference radiosonde profile in the lowest few kilometres above the surface. Sources of uncertainty were explored, and future work was suggested, such as improvements in the instrument design to minimise multipath contamination and the development of a two-dimensional adjoint model to capture horizontal variations in refractivity.

Appendix A

A more general constrained optimisation problem can be expressed in terms of the state variables v and the control variables p as

$$\begin{array}{}\text{(A1a)}& {\displaystyle }\underset{\mathit{p}}{min}f\left(\mathit{v}\left(\mathit{p}\right)\right),\text{(A1b)}& {\displaystyle }\text{subject to}\mathit{g}\left(\mathit{v},\mathit{p}\right)=\mathrm{0},\end{array}$$

where f is the scalar functional we seek to minimise and g denotes the constraints on the optimisation (equivalent to the constraint Eqs. (11b–11e) in the original optimisation problem). In a data assimilation approach, g is sometimes referred to as the forward model and describes the physics of the system. We can convert the constrained optimisation problem described by Eq. (A1) into an unconstrained optimisation problem by introducing two vectors of slack variables known as Lagrange multipliers. The Lagrangian, ℒ, can be defined in terms of f and the constraints g as

$$\begin{array}{}\text{(A2)}& \mathcal{L}\left(\mathit{v},\mathit{p},\mathbf{\Lambda }\right)\equiv f\left(\mathit{v}\right)-<\mathbf{\Lambda },\mathit{g}\left(\mathit{v},\mathit{p}\right)>,\end{array}$$

where Λ is a vector of adjoint variables associated with the constraints g and $<\cdot ,\cdot >$ is the inner product. The minimum of the Lagrangian, ℒ, and the scalar functional, f, will be equivalent when all the variables except for p are stationary points (that is, variations of the Lagrangian with respect to v and Λ are equal to 0). The total derivative of ℒ with respect to p is written as

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\mathit{p}}}={\displaystyle \frac{\partial f}{\partial \mathit{v}}}{\displaystyle \frac{\mathrm{d}\mathit{v}}{\mathrm{d}\mathit{p}}}-{\displaystyle \frac{\mathrm{d}{\mathbf{\Lambda }}^T}{\mathrm{d}\mathit{p}}}\mathit{g}-{\mathbf{\Lambda }}^T\left({\displaystyle \frac{\partial \mathit{g}}{\partial \mathit{v}}}{\displaystyle \frac{\mathrm{d}\mathit{v}}{\mathrm{d}\mathit{p}}}+{\displaystyle \frac{\partial \mathit{g}}{\partial \mathit{p}}}\right).\end{array}$$

Because the constraints g = 0, the second term on the right-hand side of Eq. (A3) can be dropped. The expression then becomes

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\mathit{p}}}={\displaystyle \frac{\mathrm{d}\mathit{v}}{\mathrm{d}\mathit{p}}}\left({\displaystyle \frac{\partial f}{\partial \mathit{v}}}-{\mathbf{\Lambda }}^T{\displaystyle \frac{\partial \mathit{g}}{\partial \mathit{v}}}\right)-{\mathbf{\Lambda }}^T{\displaystyle \frac{\partial \mathit{g}}{\partial \mathit{p}}}.\end{array}$$

In a data assimilation scheme, the Jacobian matrix $\mathrm{d}\mathit{v}/\mathrm{d}\mathit{p}$ describes the sensitivity of the trajectory of the state variables v in the forward model to changes in the control variables p. For large numbers of control variables, the cost of computing $\mathrm{d}\mathit{v}/\mathrm{d}\mathit{p}$ quickly becomes very high. Defining the adjoint variable Λ such that

$$\begin{array}{}\text{(A5)}& \mathbf{\Lambda }{\left({\displaystyle \frac{\partial \mathit{g}}{\partial \mathit{v}}}\right)}^T={\left({\displaystyle \frac{\partial f}{\partial \mathit{v}}}\right)}^T,\end{array}$$

the first term on the right-hand side of Eq. (A4) is 0, and the need to calculate the Jacobian matrix $\mathrm{d}\mathit{v}/\mathrm{d}\mathit{p}$ is avoided (Zou et al., 1997; Gill et al., 1981).

The penalty function for a single ray, j, can be augmented with the known dynamics of the system using Lagrange multipliers:

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}\mathcal{L}& =j-\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\psi }\left(\dot{h}-u\right)\mathrm{d}r-\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\mu }\\ & \left(\dot{u}-\left(\mathrm{1}-u^{\mathrm{2}}\right)\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}+{\displaystyle \frac{\mathrm{1}}{a+h}}\right)\right)\mathrm{d}r,\end{array}\end{array}$$

where $\dot{\mathit{\psi }}$ and $\dot{\mathit{\mu }}$ are the derivatives of the adjoint variables ψ and μ with respect to r, respectively. The stationary points of the Lagrangian with respect to h and u are found as

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}\mathit{\delta }\mathcal{L}{\mathrm{|}}_h& =\mathit{\delta }j{\mathrm{|}}_h-\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\psi }\mathit{\delta }\dot{h}\mathrm{d}r-\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\mu }\left(\mathrm{1}-u^{\mathrm{2}}\right)\\ & \left({\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}\right)}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}n}{\mathrm{d}{h}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{(a+h{)}^{\mathrm{2}}}}\right)\mathit{\delta }h\mathrm{d}r=\mathrm{0},\end{array}\end{array}$$

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}\mathit{\delta }\mathcal{L}{\mathrm{|}}_u& =\mathit{\delta }j{\mathrm{|}}_u+\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\psi }\mathit{\delta }u\mathrm{d}r-\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\mu }\mathit{\delta }\dot{u}\mathrm{d}r-\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathrm{2}\mathit{\mu }u\\ & \left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}+{\displaystyle \frac{\mathrm{1}}{a+h}}\right)\mathit{\delta }u\mathrm{d}r=\mathrm{0}.\end{array}\end{array}$$

Integrating Eqs. (A7) and (A8) by parts gives

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}\mathit{\delta }\mathcal{L}{\mathrm{|}}_h& =\mathit{\delta }j{\mathrm{|}}_h-{\left[\mathit{\psi }\mathit{\delta }h\right]}_{\mathrm{0}}^{r_{\mathrm{f}}}+\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\dot{\mathit{\psi }}\mathit{\delta }h\mathrm{d}r-\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\mu }\left(\mathrm{1}-u^{\mathrm{2}}\right)\\ & \left({\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}\right)}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}n}{\mathrm{d}{h}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{(a+h{)}^{\mathrm{2}}}}\right)\mathit{\delta }h\mathrm{d}r=\mathrm{0},\end{array}\end{array}$$

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}\mathit{\delta }\mathcal{L}{\mathrm{|}}_u& =\mathit{\delta }j{\mathrm{|}}_u+\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathit{\psi }\mathit{\delta }u\mathrm{d}r-{\left[\mathit{\mu }\mathit{\delta }u\right]}_{\mathrm{0}}^{r_{\mathrm{f}}}+\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\dot{\mathit{\mu }}\mathit{\delta }u\mathrm{d}r\\ & -\underset{\mathrm{0}}{\overset{r_{\mathrm{f}}}{\int }}\mathrm{2}\mathit{\mu }u\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}+{\displaystyle \frac{\mathrm{1}}{a+h}}\right)\mathit{\delta }u\mathrm{d}r=\mathrm{0}.\end{array}\end{array}$$

The adjoint variables at r = 0 are ψ_r=0 = 0 and ${\mathit{\mu }}_{r=\mathrm{0}}=\mathrm{0}$. Therefore,

$$\begin{array}{}\text{(A11)}& {\displaystyle }\mathit{\psi }{\mathrm{|}}_{r=r_{\mathrm{f}}}={\displaystyle \frac{\partial j}{\partial h}},\text{(A12)}& {\displaystyle }\mathit{\mu }{\mathrm{|}}_{r=r_{\mathrm{f}}}={\displaystyle \frac{\partial j}{\partial u}},\text{(A13)}& {\displaystyle \frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}r}}=\mathit{\mu }(\mathrm{1}-u^{\mathrm{2}})\left({\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}\right)}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}n}{\mathrm{d}{h}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{(a+h{)}^{\mathrm{2}}}}\right),\text{(A14)}& {\displaystyle \frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}r}}=\mathrm{2}\mathit{\mu }u\left({\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}h}}+{\displaystyle \frac{\mathrm{1}}{a+h}}\right)-\mathit{\psi }.\end{array}$$

Appendix B

The reported LoS AoA of the aircraft at the location of the observer is computed in terms of the geodetic Cartesian coordinates of the observer. The geodetic Cartesian coordinates of the observer (X_o, Y_o, Z_o) can be calculated using the geodetic coordinates latitude (ϕ_o), longitude (λ_o) and height (h_o) of the observer as

$$\begin{array}{}\text{(B1)}& {\displaystyle }{X}_{\mathrm{o}}=\left(N\left({\mathit{\varphi }}_{\mathrm{o}}\right)+h_{\mathrm{o}}\right)\cos \left({\mathit{\varphi }}_{\mathrm{o}}\right)\cos \left({\mathit{\lambda }}_{\mathrm{o}}\right),\text{(B2)}& {\displaystyle }{Y}_{\mathrm{o}}=\left(N\left({\mathit{\varphi }}_{\mathrm{o}}\right)+h_{\mathrm{o}}\right)\cos \left({\mathit{\varphi }}_{\mathrm{o}}\right)\sin \left({\mathit{\lambda }}_{\mathrm{o}}\right),\text{(B3)}& {\displaystyle }{Z}_{\mathrm{o}}=\left(N\left({\mathit{\varphi }}_{\mathrm{o}}\right)+h_{\mathrm{o}}\right)\sin \left({\mathit{\varphi }}_{\mathrm{o}}\right),\end{array}$$

where N is the prime vertical radius, defined in terms of the geodetic latitude ϕ as

$$\begin{array}{}\text{(B4)}& N\left(\mathit{\varphi }\right)={\displaystyle \frac{a}{\sqrt{\mathrm{1}-e^{\mathrm{2}}{\sin }^{\mathrm{2}}\left(\mathit{\varphi }\right)}}},\end{array}$$

and e² = $\mathrm{1}-b^{\mathrm{2}}/a^{\mathrm{2}}$, where a = 6378.137 km and b = 6356.75231425 km (Decker, 1986; Hofmann-Wellenhof et al., 1997). The ellipsoid Earth geometry is shown in Fig. B1.

Figure B1 The ellipsoid Earth model (thick grey line). A sphere with a radius equal to the prime vertical radius is shown.

The geocentric coordinates of the aircraft ($X_p^{\prime }$, $Y_p^{\prime }$, $Z_p^{\prime }$) are given by

$$\begin{array}{}\text{(B5)}& {\displaystyle }{X}_p^{\prime }=\left(N\left({\mathit{\varphi }}_p\right)+h_p\right)\cos \left({\mathit{\varphi }}_p\right)\cos \left({\mathit{\lambda }}_p\right),\text{(B6)}& {\displaystyle }{Y}_p^{\prime }=\left(N\left({\mathit{\varphi }}_p\right)+h_p\right)\cos \left({\mathit{\varphi }}_p\right)\sin \left({\mathit{\lambda }}_p\right),\text{(B7)}& {\displaystyle }{Z}_p^{\prime }=\left(N\left({\mathit{\varphi }}_p\right)(\mathrm{1}-e^{\mathrm{2}})+h_p\right)\sin \left({\mathit{\varphi }}_{\mathrm{o}}\right),\end{array}$$

where ϕ_p, λ_p and h_p are the geodetic latitude, geodetic longitude and altitude of the aircraft above the ellipsoid, respectively. The coordinates of the aircraft in terms of the geodetic coordinates of the observer are

$$\begin{array}{}\text{(B8)}& {\displaystyle }{X}_p=X_p^{\prime },\text{(B9)}& {\displaystyle }{Y}_p=Y_p^{\prime },\text{(B10)}& {\displaystyle }{Z}_p=Z_p^{\prime }+e^{\mathrm{2}}N\left({\mathit{\varphi }}_{\mathrm{o}}\right)\sin \left({\mathit{\varphi }}_{\mathrm{o}}\right).\end{array}$$

The geometry is shown in Fig. B2. Defining the position vectors of the observer and aircraft as ${\mathit{r}}_{\mathrm{o}}=(X_{\mathrm{o}},Y_{\mathrm{o}},Z_{\mathrm{o}})$ and ${\mathit{r}}_p=(X_p,Y_p,Z_p)$, respectively, the vector between the position of the observer and the aircraft is

$$\begin{array}{}\text{(B11)}& \mathrm{\Delta }\mathit{r}={\mathit{r}}_p-{\mathit{r}}_{\mathrm{o}}.\end{array}$$

Therefore, the zenith angle, Ω, of the aircraft as seen by the observer can be written in terms of the unit vectors:

$$\begin{array}{}\text{(B12)}& \mathrm{\Omega }={\cos }^{-\mathrm{1}}\left(\mathrm{\Delta }\widehat{\mathit{r}}\cdot {\widehat{\mathit{r}}}_{\mathrm{o}}\right).\end{array}$$

The reported LoS AoA of the aircraft is therefore equal to $\mathit{\pi }/\mathrm{2}-\mathrm{\Omega }$. The mean local radius of curvature of an ellipsoid at a given point and azimuth direction is given by (Jekeli, 2006)

$$\begin{array}{}\text{(B13)}& R_{\mathrm{roc}}={\left({\displaystyle \frac{{\sin }^{\mathrm{2}}\left(\mathit{\alpha }\right)}{N}}+{\displaystyle \frac{{\cos }^{\mathrm{2}}\left(\mathit{\alpha }\right)}{M}}\right)}^{-\mathrm{1}},\end{array}$$

where M is the radius of curvature of the meridian ellipse, given by

$$\begin{array}{}\text{(B14)}& M={\displaystyle \frac{a(\mathrm{1}-e^{\mathrm{2}})}{{\left(\mathrm{1}-e^{\mathrm{2}}{\sin }^{\mathrm{2}}\left({\mathit{\varphi }}_{\mathrm{o}}\right)\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}}}.\end{array}$$

The radius of the spherical Earth used in the ray tracing model was assumed to be the local radius of curvature of the ellipsoid at the position of the receiver, R_roc ≈ 6383.57 km (assuming an azimuth α = 45° northeast).

Figure B2 The Z coordinate of the aircraft in terms of the geodetic coordinates of the observer.