Detection and localization of F-layer ionospheric irregularities with the back-propagation method along the radio occultation ray path

The back-propagation (BP) method consists of diffractive integrals computed over a trajectory path, projecting a signal to different planes. It unwinds the diffraction and multipath, resulting in minimum disturbance to the BP amplitude when the auxiliary plane coincides with the region causing the diffraction. The method has been previously applied in global navigation satellite system (GNSS) radio occultation (RO) measurements to estimate the location of ionospheric irregularities but without complementary data to validate the estimation. In this study, the BP method is applied to a Communications/Navigation Outage Forecasting System (C/NOFS) occultation event containing scintillation signatures caused by an equatorial plasma bubble (EPB), which was parameterized with the aid of collocated data and reproduced in a wave optics propagator (WOP) simulation. In addition, a few more test cases were designed to assess the BP method with regard to the size, intensity, and placement of single- and multiple-irregularity regions. The results show a location estimate accuracy following the resolution at which the method is implemented (single bubble, reference case), whereas a bias is observed in multiple-bubble scenarios. The minimum detectable disturbance level and the estimation accuracy depend on the receiver noise level and, in the case of several bubbles, on the distance between them. These remarks provide insight into the BP results for two Constellation Observing System for Meteorology Ionosphere and Climate (COSMIC) occultation events.

1 Introduction

The Huygens–Fresnel method consists of the propagation in a vacuum of a complex wave by computing a diffractive integral of the electromagnetic (EM) field over a plane to one or multiple points in space (Sommerfeld, 1964). The direct form of the line integral is extensively combined with a wave optics propagator (WOP; Knepp, 1983) in order to obtain the EM field equivalent to the global navigation satellite system's (GNSS) complex signal sampled in the low-Earth orbit (LEO) orbit after sounding the Earth's atmosphere during a radio occultation (RO) event (Bevis et al., 1992; Kursinski et al., 1997; Gorbunov and Lauritsen, 2007).

The inverse problem, from low-Earth orbit (LEO) orbit towards the GNSS satellite, has been investigated in order to disentangle the multipath and the diffraction from the received total field and to increase the resolution of the bending-angle inversion in the lower atmosphere (Gorbunov and Gurvich, 1998a, b; Dahl Mortensen, 1998). The regions with sharp gradients in refractivity, i.e. non-homogeneities, are the source of the diffraction and multipath in amplitude and phase during the forward propagation according to Huygens' principle (Sommerfeld, 1964). The inverse form of the diffractive integral, hereafter the back-propagation (BP) method, computes the projection of the complex signal to BP planes. Ideally, the estimate location will be given at the BP plane where the BP amplitude disturbance is at its global minimum level and which coincides with the region's centre. The back-propagation field is not fully comparable to the forward field, since the back projection is performed in a vacuum; i.e. the impact height on the initial plane (boundary condition) is prolonged as straight lines to each BP plane (Gorbunov and Gurvich, 1998a).

The GNSS signal also experiences multipath and diffraction during the ionospheric propagation, where plasma irregularities above ∼ 80 km altitude are responsible for rapid fluctuations in amplitude and phase, known as ionospheric scintillation (Aarons, 1982; Yeh and Liu, 1982; Wickert et al., 2004). In the E layer (∼ 90–130 km), the regions have enhanced electron density due to the concentration of metallic ions driven by wind shear, magnetic field, gravity waves, and the influx of meteors in the Earth's atmosphere, with the main occurrence in mid-latitudes and during summer and an absence around the magnetic Equator (Arras and Wickert, 2018; Resende et al., 2018; Yu et al., 2019, 2020; Carmona et al., 2022). In the F layer, the irregularity regions in low latitudes are commonly referred to as equatorial plasma bubbles (EPBs) or equatorial spread F (ESF). The phenomenon is driven by the Rayleigh–Taylor instability mechanism, with a higher occurrence during post-sunset hours (local time), where the recombination of ions at a low altitude creates a vertical gradient in the plasma density that extends upwards to the F-region. A natural flow from the less dense (low altitudes) to denser regions (high altitudes) creates depletion areas in the form of plumes (Kelley et al., 1981; Stolle et al., 2006). The higher turbulence and gradient in density on the edges of the up-flowing bubble distorts the EM wave and eventually creates a disruption in the operation of radio frequency systems (Kelly et al., 2014). The irregularities are observed in different scale sizes (Xiong et al., 2016), and the occurrence of EPBs has shown significant seasonal, solar-cycle, and geomagnetic-activity dependence (Stolle et al., 2006, 2008; Abdu et al., 2018; Kepkar et al., 2020). In high latitudes, the occurrence of irregularity regions is not restricted in time and mostly originates from particle precipitation triggered by geomagnetic activities besides the global plasma dynamics (Jiao and Morton, 2015; Cherniak and Zakharenkova, 2016).

Following the same principle as in the lower atmosphere, the location of ionospheric irregularities in the E- and F-layers has been estimated with the BP method along the RO ray path (Gorbunov et al., 2002; Sokolovskiy et al., 2002, 2014; Cherniak et al., 2019). Back propagation has been applied to real measurements, but the estimate accuracy has been primarily assessed with WOP simulation of a generic occultation event, including a single iso- or anisotropic irregularity region modelled by one or multiple phase screens (Sokolovskiy et al., 2002). The opportunity to compare occultation measurements collocated to independent techniques must be taken to evaluate further the capabilities of the BP method in RO measurements. In Carrano et al. (2011), the scintillation pattern observed in a measurement performed by the Communications/Navigation Outage Forecasting System (C/NOFS) satellite and caused by a plasma bubble was fully modelled thanks to the parameterization of the disturbance assisted by collocated data, in which the bubble–LEO orbit satellite distance was an important variable. Further, the BP method principle has also been applied to re-scale the scintillation observed in GNSS ground receiver measurements with different frequencies and to estimate the correlation between the signals (Carrano et al., 2012, 2014).

In this study, the BP method is further assessed with WOP simulations to determine its capabilities and limitations in the context of the detection and location of F-layer irregularity regions, i.e. plasma bubbles, in RO measurements. The modelling described in Carrano et al. (2011) is considered as the initial assessment scenario of a plasma bubble in the F-region along the RO ray path, and it was used as a reference to design a few more cases with different placements, sizes, fluctuation intensities, and numbers of irregularity regions. Section 2 introduces the concept of back propagation and its equations in the scenario of an occultation event. Section 3 describes the modelling of the ionosphere and plasma bubbles in WOP simulations. Additionally, it addresses the different test cases considered in our evaluation. The simulation results are discussed in Sect. 4 and are shown to support the interpretation of two Constellation Observing System for Meteorology Ionosphere and Climate (COSMIC) measurements reported in Cherniak et al. (2019). Finally, the conclusions of the study are summarized in Sect. 5.

2 Back propagation

Assuming the scenario of a GNSS-RO simulation, the last stage of a WOP simulation takes place in a region that can be approximated to the vacuum. Therefore, the projection of the total field in the LEO orbit can be computed with the following diffraction integral (Sommerfeld, 1964):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathit{u}}_o(x,y)& =\sqrt{{\displaystyle \frac{k}{\mathrm{2}\mathit{\pi }}}}\int \mathit{u}(x,y)\cos \mathit{\xi }\\ & {\displaystyle \frac{\exp \left(ik|\mathit{r}-{\mathit{r}}_o|-i\mathit{\pi }/\mathrm{4}\right)}{|\mathit{r}-{\mathit{r}}_o{|}^{\mathrm{1}/\mathrm{2}}}}\mathrm{d}S,\end{array}\end{array}$$

where u is the total field at the last phase screen (PS), k is the wavenumber, and ξ is the angle between the normal vector to the integration plane ($\widehat{N}$) and the segment

$$\begin{array}{}\text{(2)}& |\mathit{r}-{\mathit{r}}_o|={\left({\left[{\mathit{x}}_{\mathrm{s}}-{\mathit{x}}_o\right]}^{\mathrm{2}}+{\left[{\mathit{y}}_{\mathrm{s}}-{\mathit{y}}_o\right]}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

in which the subscripts s and o stand for the coordinates on the phase screen and in the LEO orbit, and dS corresponds to the integration path, i.e. dy in this particular case. Figure 1 shows the RO geometry considered in the computation of the diffraction propagation, where the origin of the coordinate system is the Earth's centre.

Figure 1 Diffraction propagation geometry. In simulations, the forward propagation is modelled with a wave optics propagator (WOP) to consider the ionospheric irregularities. From the rightmost phase screen (PS), the complex signal is integrated along S and projected to every point in the LEO orbit.

The total field sampled in the LEO orbit (boundary condition) corresponds to the superposition of a primary and a secondary field. The primary is radiated from the GNSS satellite, whereas the secondary one results from the vibration of ions as the primary field spreads through the ionosphere. The wave field is propagated through sharp gradients in electron density inside the bubble region, which creates nonhomogeneous advances in phase around the F layer (Culverwell and Healy, 2015). As a result, rapid variations in amplitude and phase will lead to interference in the total field (focusing and defocusing), i.e scintillation (Yeh and Liu, 1982). Figure 2 illustrates the interplay of focusing and defocusing yielded by the electron density gradient, represented in terms of refractive index (n), within the irregularity patch and the resultant total field in the observational plane.

Figure 2 Illustration of the wave field generated on the GNSS transmitter (primary field) and the ionospheric mechanisms triggered by the secondary field due to the vibration of ions, followed by the overall expansion of the wave field until the observation plane (last phase screen).

By applying the same principles as in the forward propagation, it is possible to propagate the total field sampled along the LEO orbit towards the GNSS satellite. The diffraction integral for the propagation in the opposite direction, known as the back-propagation method, is written as in Gorbunov et al. (1996), Gorbunov and Gurvich (1998b), Dahl Mortensen (1998), and Sokolovskiy et al. (2002):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\mathit{u}}_b(x,y)& =\sqrt{{\displaystyle \frac{k}{\mathrm{2}\mathit{\pi }}}}\int {\mathit{u}}_o(x,y)\cos \mathit{\xi }\\ & {\displaystyle \frac{\exp \left(-ik|{\mathit{r}}_o-{\mathit{r}}_b|+i\mathit{\pi }/\mathrm{4}\right)}{|{\mathit{r}}_o-{\mathit{r}}_b{|}^{\mathrm{1}/\mathrm{2}}}}\mathrm{d}S.\end{array}\end{array}$$

However, the integration path dS is not a vertical plane as in the forward direction but rather dS=R_LEO dΩ, and the angle ξ is given between the segment $|{\mathit{r}}_o-{\mathit{r}}_{\mathrm{b}}|$ and the normal vector to the LEO orbit. The normal vector along the curved path is defined as

$$\begin{array}{}\text{(4)}& \widehat{\mathit{N}}=-\cos \mathrm{\Omega }\widehat{x}-\sin \mathrm{\Omega }\widehat{y}.\end{array}$$

Thus, the final expression for Eq. (3) is

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathit{u}}_b(x,y)& =\sqrt{{\displaystyle \frac{k}{\mathrm{2}\mathit{\pi }}}}\int {\mathit{u}}_o(x,y)\cos \mathit{\xi }\\ & {\displaystyle \frac{\exp \left(-ik|{\mathit{r}}_o-{\mathit{r}}_b|+i\mathit{\pi }/\mathrm{4}\right)}{|{\mathit{r}}_o-{\mathit{r}}_b{|}^{\mathrm{1}/\mathrm{2}}}}{R}_{\mathrm{LEO}}\mathrm{d}\mathrm{\Omega },\end{array}\end{array}$$

where $\cos \mathit{\xi }=\widehat{N}\cdot \widehat{\mathit{r}}$. Figure 3 shows the geometry of the back-propagation scenario, the relation between the angles, and the normal vector direction that changes along the LEO orbit.

Figure 3 Back-propagation geometry. The diffraction integral is computed along the LEO orbit path (S) and projected to every point on the closest PS. Once the back-propagated signal is available on the auxiliary plane, the Zverev transform is computed iteratively towards the GPS satellite.

Slightly different procedures are described in the literature to obtain the total field at different BP planes as the complex signal is propagated towards the GNSS satellite. They are outlined as follows:

1. Compute Eq. (5) along the LEO orbit to obtain the BP signal at different BP planes (Sokolovskiy et al., 2002).

2. Compute the Zverev transform to obtain the BP signal at different BP planes (Gorbunov et al., 2002; Gorbunov and Lauritsen, 2007), which consists of propagating the BP signal at the right-most PS (via diffraction integral) by applying direct and inverse Fourier transforms recursively, viz.

$$\begin{array}{}\text{(6)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{u}}}_{\mathrm{b}}\left(x_{\mathrm{b}},k_y\right)=\mathfrak{F}\left\{{\mathit{u}}_{\mathrm{b}}\left(x_{\mathrm{b}},y_{\mathrm{b}}\right)\right\},\text{(7)}& {\displaystyle }{\mathit{u}}_{\mathrm{b}}(x,y)={\mathfrak{F}}^{-\mathrm{1}}\left\{{\stackrel{\mathrm{̃}}{\mathit{u}}}_{\mathrm{b}}\left(x_{\mathrm{b}},k_y\right)\exp \left(i\sqrt{k^{\mathrm{2}}-{\mathit{k}}_y^{\mathrm{2}}}\left(x_o-x_{\mathrm{b}}\right)\right)\right\},\end{array}$$

where 𝔉 is the Fourier operator, and k_y is the spatial angular frequency.

In order to locate the source of the secondary field, i.e. irregularity regions, the auxiliary plane with minimum BP amplitude disturbance should be found, as the disturbance is expected to reduce gradually up to the source point. The standard deviation of the BP amplitude is the metric used to quantify the disturbance in each auxiliary plane, and it is defined as in Gorbunov et al. (1996):

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_u=\sqrt{{\displaystyle \frac{{\left(\mathit{u}{}^{\prime }-\stackrel{\mathrm{‾}}{u{}^{\prime }}\right)}^{\mathrm{2}}}{N}}},\end{array}$$

where $\mathit{u}{}^{\prime }=\mathit{u}-\stackrel{\mathrm{‾}}{u}$ corresponds to detrended BP amplitude achieved by a three-pass second-order Savitzy–Golay filter and assuming a window length of 10 km according to the typical irregularity outer scale (Carrano et al., 2011; Zeng et al., 2019).

The ability to find the origin of the secondary field along the ray path is dependent on the secondary-field amplitude (proportional to the electron density gradient) and on the noise level of the LEO receiver. These aspects are investigated in simulations, including the modelling of ionospheric disturbances.

3 Ionospheric simulation

The effects of ionospheric refractivity are accounted for in a WOP simulation by assuming that the electron density profile (EDP) is a part of the atmospheric model. The refractive index, combining the neutral atmosphere and the ionosphere, is defined as

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathit{n}}_i=-\mathrm{40.3}{\displaystyle \frac{\mathit{\rho }}{f^{\mathrm{2}}}},\text{(10)}& {\displaystyle }\mathit{n}={\mathit{n}}_n+{\mathit{n}}_i,\end{array}$$

where f is the carrier frequency, ρ is the electron density (el m⁻³), and subscripts n and i denote the neutral atmosphere and ionosphere, respectively. The addition of the ionospheric model includes the respective phase shift into the total phase accumulated during the forward wave propagation. From the RO perspective, the excess path due to the ionospheric propagation in such a scenario may result in an extra accumulated bending angle proportional to f⁻² – i.e. the lower the frequency, the larger the bending. Additionally, the signals in different frequencies have different bending angles due to slightly different propagation paths. Consequently, the signals have different integrated electron densities (Culverwell and Healy, 2015).

3.1 F-region irregularity: plasma bubbles

Under low ionospheric activity, EDPs tend to resemble a slow function (Culverwell and Healy, 2015). Under high-activity periods and during the transition between day- and nighttime, there is a higher incidence of regions of localized irregularities, e.g. plasma bubbles, leading to a sharper gradient in electron density (Jiao and Morton, 2015; Kepkar et al., 2020). Such regions are responsible for large-, medium-, and small-scale irregularities, corresponding to sizes up to the Fresnel scale (Xiong et al., 2016). In an RO geometry and especially in the range of ionospheric altitudes where the bending is significantly smaller than in neutral atmosphere (Kursinski et al., 1997), the Fresnel scale is given by

$$\begin{array}{}\text{(11)}& {\displaystyle }{d}_F=\mathrm{2}\sqrt{{\displaystyle \frac{\mathit{\lambda }{L}_t{L}_r}{L_t+L_r}}},\text{(12)}& {\displaystyle }{d}_F\approx \mathrm{1.5}\mathrm{km},\end{array}$$

where λ is L1 band wavelength, L_t is the horizontal distance between the GNSS satellite and the Earth's limb ($\sim \mathrm{28.5}\times {\mathrm{10}}^{\mathrm{6}}$ m), and L_r is the LEO horizontal distance ($\sim \mathrm{3.4}\times {\mathrm{10}}^{\mathrm{6}}$ m, assuming an altitude orbit of 820 km). The propagation through these irregularities results in diffraction and refraction of the electromagnetic field. These effects are observed as abrupt fluctuations in amplitude and phase, referred to as scintillations (Aarons, 1982; Yeh and Liu, 1982; Wickert et al., 2004; Zeng and Sokolovskiy, 2010). Moreover, the presence of plasma bubbles introduces asymmetries between the inbound (GNSS to tangent point) and outbound (tangent point to LEO) segments of the ray trajectories. This condition contradicts the assumption of the spherical symmetry of the atmosphere in retrievals via the Abel transform (Fjeldbo et al., 1971), and it is related to high-order terms composing the bias after the standard ionospheric correction (Vorob'ev and Krasil'nikova, 1994). The high-order bias, critical in meteorological and climate applications, is handled by either kappa or bi-local correction (Healy and Culverwell, 2015; Liu et al., 2020).

3.1.1 Single bubble

The location estimation of the plasma bubbles in the F region is a complicated task in RO measurements. The ray path between GNSS and LEO satellites includes ionospheric propagation in two segments, i.e. ray inbound and outbound. The disturbance that was observed in the sampled signal and that originated during either the former or the latter segment cannot be visually distinguished. The back-propagation (BP) method has been used to detect irregularities in the F region in studies using both simulations and real occultation measurements (Sokolovskiy et al., 2002, 2014; Cherniak et al., 2019).

However, there is a lack of RO events combined with collocated data provided by different systems where the true location of the irregularity region is precisely known. In this study, the model of isotropic irregularities representing a plasma bubble in the equatorial region is considered in WOP simulations to evaluate the estimation obtained with the BP method. The model has been described in Carrano et al. (2011) and corresponds to a measurement performed by the C/NOFS satellite and collocated with an incoherent scatter radar and a very-high-frequency (VHF) ground-based receiver. The collocated data allowed a good estimation of the placement and size of the bubble to be made; this is in addition to the parameters required in the modelling of the disturbance observed in the occultation measurement.

The plasma bubble is modelled by a 2-D random realization of Gaussian variables filtered by the spectral density function (SDF):

$$\begin{array}{}\text{(13)}& {\mathbf{\Phi }}_{\mathrm{\Delta }\mathit{\rho }}=\mathrm{4}\mathit{\pi }{k}_{\mathrm{0}}^{\left(\mathrm{2}\mathit{\nu }-\mathrm{2}\right)}{\displaystyle \frac{\mathrm{\Gamma }\left(\mathit{\nu }\right)}{\mathrm{\Gamma }(\mathit{\nu }-\mathrm{1})}}{\displaystyle \frac{\mathrm{1}}{{\left(k_{\mathrm{0}}^{\mathrm{2}}+{\mathit{k}}_x^{\mathrm{2}}+{\mathit{k}}_y^{\mathrm{2}}\right)}^{\mathit{\nu }}}},\end{array}$$

where k_x,y are the wave numbers along and across the propagation direction, $k_{\mathrm{0}}=\mathrm{2}\mathit{\pi }/L_{\mathrm{0}}$ is the outer-scale wave number, Γ is the Euler's gamma function, and ν denotes the spectral slope. The filtered variables,

$$\begin{array}{}\text{(14)}& \mathbf{\Delta }\mathit{\rho }={\mathfrak{F}}^{-\mathfrak{1}}\left\{\sqrt{{\mathbf{\Phi }}_{\mathrm{\Delta }\mathit{\rho }}\mathrm{SF}}{\mathbf{r}}_{\mathrm{m}}\right\},\end{array}$$

are modulated to the electron density model,

$$\begin{array}{}\text{(15)}& \mathit{\rho }={\mathit{\rho }}_{\mathrm{b}}\left[\mathrm{1}+\mathbf{\Delta }\mathit{\rho }{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}\mathit{B}\right],\end{array}$$

where ρ_b is the background EDP, and ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}$ is the root mean square (rms) level of the fluctuations. The bubble width is a Gaussian envelope function,

$$\begin{array}{}\text{(16)}& {\displaystyle }\mathit{B}(x,y)=e^{\frac{{\left[\mathit{\alpha }-{\mathit{\alpha }}_{\mathrm{0}}\right]}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}}},\text{(17)}& {\displaystyle }\mathit{\alpha }(x,y)={\tan }^{-\mathrm{1}}[\mathit{y}/(x+R_{\mathrm{e}}\left)\right],\end{array}$$

in which the function maximum and the bell width are set by

$$\begin{array}{}\text{(18)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{0}}={\tan }^{-\mathrm{1}}\left({\displaystyle \frac{x_{\mathrm{0}}}{hmF\mathrm{2}+R_{\mathrm{e}}}}\right),\text{(19)}& {\displaystyle }{\mathit{\sigma }}_{\mathit{\alpha }}={\displaystyle \frac{L_{\mathrm{H}}}{A(hmF\mathrm{2}+R_{\mathrm{e}})}},\end{array}$$

where x₀ denotes the bubble placement along the x axis, hmF2 is the F-region electron density peak height, L_H corresponds to the bubble width, R_e is the Earth's radius, and the scaling factor A≈1.348. In Eq. (14), r_m corresponds to the grid of Gaussian random numbers, and $\mathrm{SF}=L/\mathrm{2}\mathit{\pi }$ is a spatial factor in which L is the bubble vertical extension.

The set of parameters estimated in Carrano et al. (2011) was used in our WOP simulation to replicate the scintillation in the total field with equivalent deterministic properties. Further details about the implementation of the wave optics propagator used in the simulations are given in Ludwig-Barbosa et al. (2020) and Ludwig‐Barbosa et al. (2020). Figure 4 shows the Gaussian envelope and the filtered random realization modulated to the electron density model.

Figure 4 Modelling of the irregularity region used during the forward propagation in simulations. (a) Gaussian envelope defined by Eq. (16). (b) Resultant irregularities modulated to a background EDP, as given in Eq. (15).

Figure 5 shows the normalized signal intensity at the observational plane and the power spectral density (PSD) computed within 280 and 340 km. The results are in good agreement with the ones reported in Carrano et al. (2011) and validate our WOP simulation.

Figure 5 WOP results considering the set of parameters described in Carrano et al. (2011). The original C/NOFS measurement had an average SNR level of around 1500 V V⁻¹.

The simulated total field disturbed by the plasma bubble during the forward propagation, and sampled at the right-most PS (see Fig. 1), is considered to be the boundary condition of the BP method. The scenario is used as the reference model for different test cases to assess the capabilities and limitations of the BP method in the presence of a single plasma bubble, namely

• The accuracy of the location estimate along the x axis. the position of the region of irregularities is controlled by modifying x₀ in Eq. (16).

• The accuracy of the location estimate along the x axis with different rms fluctuation levels. The level of irregularities modulated to the EDP is defined by ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}$ in Eq. (15).

• The accuracy of the location estimate along the x axis with different vertical extensions of the bubble.

• The accuracy of the location estimate along the x axis with different bubble width. The extension along x axis is controlled by L_H in Eq. (19).

3.1.2 Multiple bubbles

In addition to the single-bubble cases, a second plasma bubble was added to the ray trajectory by superposing another envelope function on the one shown in Fig. 4a by simply assuming a different x₀ in Eq. (16). The test cases with two plasma bubbles allow us to evaluate the BP method in the following scenarios:

• the accuracy of the location estimate along the x axis for the two plasma bubbles

• the accuracy of the location estimate in relation to different separation distances between bubbles

• the accuracy of the location estimate when bubbles have different rms fluctuation levels.

Figure 6 Scenario with a single F-layer bubble in the inbound region. The irregularity region centre is around $x_{\mathrm{ref}}=-\mathrm{346.7}$ km, assuming ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$. Vertical black lines represent the detrended BP amplitudes computed in each auxiliary plane (50 km resolution representation). The blue curve corresponds to the detrended BP amplitude standard deviation, computed every 5 km. The minimum σ_u provides the estimated centre of the irregularity patch. Estimation error ${\mathit{\epsilon }}_x=-\mathrm{1.7}$ km.

4 Results

In the simulations, the filtered random field was modulated with an EDP modelled by Chapman's function (Culverwell and Healy, 2015) considering the F-region peak ($nmF\mathrm{2}=\mathrm{8.81}\times {\mathrm{10}}^{\mathrm{11}}$ el m⁻³), height (hmF2=288.5 km), and scale height (H=31 km) according to the EDP described in Carrano et al. (2011).

The forward-propagation simulations of the test cases did not include the propagation to the LEO orbit via the diffraction integral in Eq. (1). Therefore, the BP signals are computed via Eqs. (6), (7) since the boundary condition is given on the vertical plane. The WOP signals include instrument noise, which assumed a Meteorological Operational satellite (MetOp-A) occultation event at a low latitude as a reference to the signal-to-noise (SNR) level (see Fig. A1 in the Appendix). This measurement extends up to 600 km straight-line tangent altitude (SLTA), an exceptional feature compared to nominal MetOp measurements. Normally, the GNSS signal is tracked up to around 100 km SLTA, but an experimental campaign during the MetOp-A end-of-life operation had its tracked SLTA range extended to the point where the F region was included. Differently from the neutral atmosphere region, the SNR level decays with altitude due to the antenna radiation pattern. At this particular measurement (and in simulations), the SNR reference level assumed to estimate the instrumental noise strength in the F-region peak was around 600 V V⁻¹.

In forward-propagation simulations, the closest phase screen to x₀, defining the centre of the irregularity region, was placed at −346.7 km. Therefore, this was the placement reference (x_ref) assumed in the accuracy analysis. The BP planes were computed every 5 km, which defines the precision of the estimations in our implementation.

Figure 7 Scenario with single F-layer bubble in the outbound region. The irregularity region was centred around x_ref=346.7 km, assuming ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$. Estimation error ε_x=1.7 km.

4.1 Single bubble

Figure 6 shows the BP amplitudes when bubbles were placed at $x_{\mathrm{ref}}=-\mathrm{346.7}$ km, with the assumed rms fluctuation level ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$ (Carrano et al., 2011). The plasma bubble is represented in the background of the BP amplitudes (black vertical lines).

The rms fluctuation level corresponds to a variation of ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}\approx \pm \mathrm{1.5}\times {\mathrm{10}}^{\mathrm{11}}{\mathrm{el}\mathrm{m}}^{-\mathrm{3}}$, which results in weak scattering (S₄=0.26) in agreement with Carrano et al. (2011). The estimate error corresponds – i.e. ${\mathit{\epsilon }}_x=x_{\mathrm{ref}}-x_{{\mathit{\sigma }}_u,\min }=-\mathrm{1.7}$ km.

Figure 7 shows the result considering the bubble placement on the outbound ray path.

In the single-bubble scenario, the location estimate has good accuracy regardless of the placement in the inbound or outbound sector. Therefore, the location estimate in a single-bubble scenario is limited by the precision considered in the BP method, herein 5 km. The minor difference in the scintillation index (S₄) is related to the filtered random variables assumed in the bubble model, which can create a variation in the resultant electron density obtained in the simulation.

4.1.1 Influence of rms fluctuation level

A parametric evaluation of ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}$ was performed to assess the minimum fluctuation level at which the bubble is detectable with the BP method. Figure 8 shows the box chart comparing the sensitivity in detection for the three different levels, ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{2}$ %, 3.0 %, and 17 % (reference case) in terms of estimation accuracy along the x axis, and the correspondent BP standard deviation curves.

Figure 8 Influence of fluctuation level in simulations assuming a single irregularity region in the inbound region. (a) Box chart comparing different rms fluctuation levels. The vertical dashed line indicates the bubble placement. The detection improves significantly when the fluctuations level is ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}>\mathrm{3.0}$ % ($\pm \mathrm{2.64}\times {\mathrm{10}}^{\mathrm{10}}$ el m⁻³), with an estimate median $\stackrel{\mathrm{‾}}{x}=-\mathrm{327.5}$ km achieved with such a set-up. Weaker irregularities, e.g. ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{2}$ %, are not distinguishable from the receiver noise and yield a poor location estimate of the irregularity patch. (b) BP amplitude standard deviation. Shaded regions depict the 2σ interval. The same colour scheme is used in both figures, and the grey line represents the receiver noise level. Figures depict results assuming 20 realizations for each fluctuation level.

The curve corresponding to ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}\le \mathrm{2}\mathit{\%}$ (red curve in Fig. 8b) does not have a clear global minimum. The BP standard deviation level lies beneath the threshold value determined by the receiver noise level around hmF2 (σ₀≈0.0456 after Figs. 8b and A1b). Thus, this indicates that the estimations are unreliable when σ_u≤σ₀. For rms fluctuation levels ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}>\mathrm{3.0}\mathit{\%}$ ($\pm \mathrm{2.64}\times {\mathrm{10}}^{\mathrm{10}}$ el m⁻³), the region of irregularities is detectable with a median of $x=-\mathrm{327.5}$ km and the interval of $[-\mathrm{410},-\mathrm{200}]$ km corresponding to 50 % of the estimates ($\mathrm{63.3}<{\mathit{\epsilon }}_x<-\mathrm{146.7}$ km). Regardless, the disturbance level ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{3.0}\mathit{\%}$ assumed in the forward propagation corresponds to a weal disturbance in the LEO's orbit (S₄<0.1). At this level of scintillation, the disturbance created by ionospheric irregularities cannot be distinguished from other sources of error. Therefore, the low accuracy achieved in the estimation is not a concerning result (Béniguel et al., 2009; Ma et al., 2019). For the reference case, an accuracy around the method precision was achieved $({\mathit{\epsilon }}_x=-\mathrm{3.3}$ km).

4.1.2 Influence of bubble vertical extension

Figure 9 shows the comparison between the location estimate obtained with WOP simulations assuming different vertical thicknesses for the irregularity region and ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$. The thickness was controlled by applying a Tukey window to the right term in Eq. (15).

The dashed black curve shows the standard deviation curve for the bubble with original dimensions, in which the effective extension of the bubbles is defined by the region around the F region, with an electron density within 75 % of the peak value (∼ 60 km; Carrano et al., 2011). The maximum estimation error observed in simulations was −1.7 km for all cases. This result implies that the vertical extension of the region does not impact the location estimate.

Despite the disturbance in the simulation being located in the F region, the vertical extensions shorter than 10 km resemble the thickness of a sporadic E layer (Zeng and Sokolovskiy, 2010; Arras and Wickert, 2018), and it confirms the capability presented in Gorbunov et al. (2002), Sokolovskiy et al. (2014), and Cherniak et al. (2019). Moreover, the scintillation in the E layer may have a potential advantage in the purview of accuracy given the higher SNR level (lower noise floor) around 100 km (see Fig. A1b).

Figure 9 BP amplitude standard deviation in scenarios assuming different vertical extensions of a single bubble in the inbound region and assuming ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$. The dashed black curve corresponds to the reference case. The location estimate is possible up to the thinnest layer, resembling the sporadic E-layer dimension. Estimation error ${\mathit{\sigma }}_x=-\mathrm{1.7}$ km. Vertical dashed lines indicate the placement of the irregularity patch.

4.1.3 Influence of bubble width

Figure 10 shows results for scenarios assuming different bubble widths and fixed fluctuation levels (${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$). A region with extension L_H≤20 km creates low scintillation in the GNSS signal (S₄<0.2), but it is still detectable, and it has estimation error ${\mathit{\epsilon }}_x=-\mathrm{5}$ km. Narrower regions do not show a clear global minimum, since the disturbances are at the same level as the receiver noise.

Figure 10 Single bubble with different widths (L_H). The detection is possible for wide regions, but estimate accuracy decreases with increasing width.

The detection of irregularities is theoretically possible even for wide regions, which leads to higher disturbances, as indicated by the scintillation index. However, the uncertainty about its centre estimate increases proportionally with the region width despite the increasing difference between the global minimum level and the noise floor. Thus, the extension of the irregularity region must be shorter than the distance between GNSS and LEO satellites, as stated in Sokolovskiy et al. (2002).

Figure 11 Bubbles at $x_{\mathrm{ref},\mathrm{1}}=-\mathrm{346.7}$ and $x_{\mathrm{ref},\mathrm{2}}=\mathrm{346.7}$ km; ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$. Estimation error ε_x2=71.7 km.

4.2 Multiple bubbles

Figure 11 shows two bubbles symmetrically placed around the origin at $x_{\mathrm{ref},\mathrm{1}}=-\mathrm{346.7}$ and $x_{\mathrm{ref},\mathrm{2}}=\mathrm{346.7}$ km and with the same fluctuation level (${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$).

The global minimum σ_u corresponds to the bubble placed on the outbound region, the last irregularity region along the ray path (forward propagation). The accuracy of the location estimate is affected significantly by the presence of the inbound bubble and by the instrument noise, yielding an estimation error ε_x≈71.7 km. The location estimate of the inbound bubble is a rather complicated task, since the presence of the outbound bubble shadows its contribution to the total wave field; therefore, a clear local minimum is not detectable in the BP amplitude standard deviation.

Figure 12 shows the scenario with a larger separation between the irregularity regions (Δx=1200 km). The minima become more distinguishable and improve their location estimates slightly. Nevertheless, the most accurate estimation is given on the outbound bubble (ε_x≈40 km), with the instrument noise having a partial contribution to the error. Regarding the inbound bubble, there is an indication of the irregularity placement around $x=-\mathrm{500}$ km, but the estimation error is greater than for the outbound patch.

Figure 12 Bubbles at $x_{\mathrm{1}}=-\mathrm{600}$ and x₂=600 km; ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$. Estimation error ε_x2=40 km. Note the indication of the inbound bubble's placement along the ray path.

A comparison between Figs. 10c, d (wide-bubble scenarios) and Figs. 11 and 12 shows that it is possible to distinguish cases with a single wide irregularity region from a scenario with multiple smaller bubbles, since the latter would likely present more than one local minimum along the ray path. Nevertheless, the location estimates of secondary patches are less reliable.

In contrast to single-region cases where the predominant constraint to detection is the noise level (σ_u≈σ₀), these results indicate that the separation between the regions has a major influence on the detection and location task of multiple patches.

Figure 13 Bubbles at $x_{\mathrm{ref}}=\pm \mathrm{346.7}$ km with weaker rms fluctuation levels on the inbound (a, b) and the outbound bubble (c, d).

4.2.1 Influence of rms fluctuation level

In these test cases, the rms fluctuation level of one of the bubbles was kept constant (${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$), while the other had the fluctuation set to weaker values. Figure 13 depicts the results assuming ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{6}\mathit{\%}$ and ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{12}\mathit{\%}$.

The standard deviation curves in scenarios including a bubble with ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{6}\mathit{\%}$ are similar to the one observed in the scenario of a single bubble (see Figs. 6 and 7). However, the location of the global minima along the x axis differs, indicating that the presence of a weaker bubble affects the location estimate of the predominant irregularity region. The remarks are valid despite the placement of the weaker region in the inbound or outbound sector. However, the inbound bubbles have a greater impact on the estimation accuracy of the inbound bubbles than in the opposite case.

Figure 14 shows the comparison of standard deviation curves assuming different rms fluctuation levels on the bubble placed at the inbound sector.

Figure 14 Comparison between simulations assuming different rms fluctuation levels on the inbound bubble and constant fluctuation levels on the outbound bubble (${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho },\mathrm{2}}=\mathrm{17}\mathit{\%}$). The legend shows the ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho },\mathrm{1}}$ level and the location estimate along x, given after the σ_u global minimum. Vertical dashed lines depict the placement of the irregularity regions in the simulations. The dashed black curve corresponds to the case of a single bubble in the outbound sector. The dashed blue curve corresponds to the scenario shown in Fig. 11.

Figure 15 BP amplitudes of COSMIC occultation events during geomagnetic storms. The right y axis corresponds to the values of the normalized BP amplitude standard deviation (red curve). The global minimum in the red curves estimates a similar position of the main irregularity region for both occultations (x=2600–2800 km). The arrows highlight local minima, which may indicate the existence of other regions of ionospheric irregularities. Measurement average SNR levels of (a) 564 V V⁻¹ and (b) 694 V V⁻¹.

A clear shift of the global minimum towards the weaker patch is observed around x_ref as ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho },\mathrm{1}}$ increases from 6 % to 17 % (${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho },\mathrm{1}}={\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho },\mathrm{2}}$), which leads to a gradual increase in the estimation error. After ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho },\mathrm{1}}>{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho },\mathrm{2}}$, the estimation indicates the position of the inbound bubble and no longer the outbound one.

4.3 Analysis of COSMIC occultation events

The remarks made after the test cases are used in the evaluation of two COSMIC occultation events presented in Cherniak et al. (2019). The measurements were performed during a severe geomagnetic storm between 22–23 June 2015. Their results are replicated in Fig. 15 after using Eq. (5) to compute the BP amplitude at x=3000 km, followed by employing Eqs. (6), (7) recursively to obtain the total field at the other auxiliary planes.

The global minima are found between 2600–2800 km in both occultations, indicating the position of the main region of irregularities along the ray path. In Fig. 15a, the BP amplitude standard deviation was computed assuming the entire height range available in every BP plane, since the measurement SNR (figure not shown) did not contain any clear signature of sporadic-E scintillation (Zeng and Sokolovskiy, 2010; Arras and Wickert, 2018; Yu et al., 2020; Carmona et al., 2022). In Fig. 15b, a u-shaped fade was presented around 100 km SLTA (figure not shown). This altitude corresponds to the conventional range of occurrence for sporadic Es, likely indicating that the irregularities were aligned with the propagation direction. Therefore, the height range around the u-shaped fade has not been included in the calculation of the BP amplitude standard deviation, and so it does not affect the location estimate of disturbances in the F region. The same methodology has been assumed in Gorbunov et al. (2002) and Cherniak et al. (2019).

As seen in simulations, the existence of other local minima in the standard deviation curve gives the indication of not only one but two or more irregularity regions during the occultation events. The arrows point out the approximate location of these local minima in Fig. 15. The confirmation of the existence of such regions requires collocated measurements, similar to the case reported in Carrano et al. (2011). In addition, the existence of multiple regions has been shown to reduce the estimate accuracy given after the global minimum to some extent. In reality, the main irregularity region could have been placed slightly closer to the LEO satellite than the position estimated by the BP method, whereas the secondary patches may be slightly farther away from the receiver (similar to Fig. 12).

5 Conclusions

The capability of back propagation to detect irregularity regions in the F layer, e.g. ionospheric plasma bubbles, has been assessed with WOP simulations. The reference case corresponded to a single bubble at the inbound sector observed in a C/NOFS occultation event, in which the location, size, and distance from the LEO orbit have been confirmed with collocated data (Carrano et al., 2011). The same model of isotropic irregularities was applied to all the other test scenarios evaluated with WOP simulations.

In the simulation of single-bubble scenarios, the location estimate accuracy of the irregularity region along the ray path follows the method resolution for the reference case (${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\rho }/\mathit{\rho }}=\mathrm{17}\mathit{\%}$). The bubble placement in either inbound or outbound regions did not affect the detection and location estimate of the irregularity regions. Additionally, the detection of bubbles has been possible regardless of the region width or vertical extension when S₄>0.2. However, the accuracy of the centre estimate decreases with increasing width.

In multiple-bubble scenarios, the ability to estimate the location of bubbles requires the patches to be well separated. Then, the regions are detectable, but the accuracy of the estimate differs. The region yielding the stronger disturbance (predominant) has the most accurate location estimate. However, a bias towards the weaker bubbles is inherent, and it increases with the rms fluctuation level. If secondary bubbles have a very weak fluctuation strength, the patch is shadowed by the dominant region, and their existence can be untraceable. In the case of aligned bubbles with similar intensities, the most accurate estimation corresponds to the latest region along the forward-propagation direction.

Most importantly, the capability of the detection and location of irregularity patches has been shown to be limited by the receiver noise level; i.e. localizing irregularity patches with the BP method is unfeasible when the noise level is greater than the amplitude of the ionospheric scintillation (σ_u<σ₀). At the SNR level assumed as the reference in our simulations (MetOp), even irregularity patches in the F region corresponding to low scintillation were detectable ($\mathrm{\Delta }\mathit{\rho }\sim \pm \mathrm{2.64}\times {\mathrm{10}}^{\mathrm{10}}$ el m⁻³). This fluctuation corresponds to the local gradient within the bubble region and, therefore, depends on the local mean density (background EDP), the patch size, and the distance between the bubble and receiver. Nevertheless, the minimum detectable level will vary among different receivers according to their noise figures.

The SNR levels, as well as the highest SLTA points in measurements, differ in different RO missions. An SLTA range which includes the ionospheric layer, i.e. further than 100 km SLTA, as seen in the experimental MetOp-A campaign, is an important feature to accurately detect and locate the ionospheric plasma bubbles in RO measurements. A minimized influence of the antenna gain in higher SLTA might also contribute to improving the results obtained with the BP method. Nonetheless, the results indicate that the present operating SNR level in the MetOp constellation is sufficient to detect even low scintillation levels.

The information about the location of irregularity regions, e.g. plasma bubbles, is relevant in ionospheric modelling and could potentially support mapping such phenomena and their climatology. In this context, RO data have the potential to improve the gaps in the coverage provided by networks of ground-based receivers that detect and track these regions. Our results should be taken as a complement to the investigations described in Gorbunov et al. (2002), Sokolovskiy et al. (2002), and Cherniak et al. (2019). Further evaluations of collocated occultation events with data provided by different systems, in line with Carrano et al. (2011), are required to evaluate the method capabilities (also with regard to E layers). In combination with the location along x, the horizontal and vertical extension of the plasma irregularity are also parameters of great interest for modelling the plasma irregularities. Approaches to estimating such features, as well as alternatives to locate secondary regions, should be investigated in future works.

Appendix A: Including instrument noise in WOP signals

Figure A1 (a) L1 coarse acquisition (C/A) code SNR_1Hz: the blue curve shows the original SNR, and the red curve depicts the averaged curve, the values of which were used as a reference in Eq. (A18). The decay in SNR observed with increasing SLTA (> 100 km) is due to the antenna gain pattern. (b) The amplitude of the noise added to the WOP signal. (c) WOP amplitude with and without added noise on the observational plane (last PS; single-inbound-bubble scenario) (Carrano et al., 2011).

In WOP simulations, the signal transmitted by the GNSS (boundary condition) is assumed to be a cylindrical wave. The propagation between the GNSS satellite and the first phase screen occurs in free space, with amplitude decay $\propto \mathrm{1}/\sqrt{r}$. For the sake of practicality, the complex signal is normalized on the first PS. Then, the medium refractivity is recursively accounted for by modifying the instant phase of the incident wave and propagating it in a vacuum until the next neighbouring phase screen (Knepp, 1983). At the last PS, the normalized complex signal in the WOP ($\widehat{u}$) can model a real signal by using a constant calibration factor, A, viz.

$$\begin{array}{}\text{(A1)}& u_{\mathrm{signal}}\left(t\right)=A\widehat{u}\left(t\right).\end{array}$$

The total signal will also include noise,

$$\begin{array}{}\text{(A2)}& u_{\mathrm{total}}\left(t\right)=u_{\mathrm{signal}}\left(t\right)+u_{\mathrm{noise}}\left(t\right).\end{array}$$

We used the measured SNR from a representative MetOp-A occultation event to estimate the appropriate noise level added to the WOP amplitude,

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\widehat{u}}_{\mathrm{noise}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{A}}{u}_{\mathrm{noise}}\left(t\right),\text{(A4)}& {\displaystyle }{\widehat{u}}_{\mathrm{total}}\left(t\right)=\widehat{u}\left(t\right)+{\widehat{u}}_{\mathrm{noise}}\left(t\right).\end{array}$$

The noise in occultation measurements has several sources: thermal noise in the receiver, clock noise, and co-channel noise. For this task, we assumed a normal distribution to model the white noise – i.e. $X,Y\sim \mathcal{N}(\mathit{\mu },{\mathit{\sigma }}^{\mathrm{2}})$, where μ is the mean value, and σ² is the variance. Then, the noise in the ith sample is

$$\begin{array}{}\text{(A5)}& {\displaystyle }{u}_{\mathrm{noise}}\left(t_i\right)={\mathit{\sigma }}_{\mathrm{0}}(X+jY)/\sqrt{\mathrm{2}},\text{(A6)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{0}}{}^{\prime }={\mathit{\sigma }}_{\mathrm{0}}/A,\text{(A7)}& {\displaystyle }{\widehat{u}}_{\mathrm{noise}}\left(t_i\right)={\mathit{\sigma }}_{\mathrm{0}}{}^{\prime }(X+jY)/\sqrt{\mathrm{2}},\end{array}$$

where $X,Y\sim \mathcal{N}(\mathrm{0},\mathrm{1})$. Next, we obtain the average noise power (approximation due to the finite number of samples) by multiplying the noise with its complex conjugate and taking the average over a large time window,

$$\begin{array}{}\text{(A8)}& P_{\mathrm{noise}}=〈u_{\mathrm{noise}}{u}_{\mathrm{noise}}^{\ast }〉\approx {\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}.\end{array}$$

Likewise, the averaged signal power becomes

$$\begin{array}{}\text{(A9)}& P_{\mathrm{signal}}=〈u_{\mathrm{signal}}{u}_{\mathrm{signal}}^{\ast }〉\approx A^{\mathrm{2}}〈\widehat{u}{\widehat{u}}^{\ast }〉.\end{array}$$

The SNR in terms of the signal and noise power, with units [W W⁻¹], is given by

$$\begin{array}{}\text{(A10)}& {\mathrm{SNR}}_{\mathrm{W}}={\displaystyle \frac{P_{\mathrm{signal}}}{P_{\mathrm{noise}}}}\approx {\displaystyle \frac{A^{\mathrm{2}}〈\widehat{u}{\widehat{u}}^{\ast }〉}{{\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}}}.\end{array}$$

Hence,

$$\begin{array}{}\text{(A11)}& {\mathit{\sigma }}_{\mathrm{0}}=\sqrt{{\displaystyle \frac{A^{\mathrm{2}}〈\widehat{u}{\widehat{u}}^{\ast }〉}{{\mathrm{SNR}}_{\mathrm{W}}}}},\end{array}$$

and

$$\begin{array}{}\text{(A12)}& {\mathit{\sigma }}_{\mathrm{0}}{}^{\prime }=\sqrt{{\displaystyle \frac{〈\widehat{u}{\widehat{u}}^{\ast }〉}{{\mathrm{SNR}}_{\mathrm{W}}}}}.\end{array}$$

In case different sample rates are used in the measurements and the simulations, one has to take into account the sample rate or the bandwidth (B), where

$$\begin{array}{}\text{(A13)}& B\propto f_{\mathrm{s}},\end{array}$$

in which f_s is the sample rate in Hz. The noise power is given by

$$\begin{array}{}\text{(A14)}& P_{\mathrm{noise}}=B{N}_{\mathrm{0}},\end{array}$$

where N₀ is the noise power density in W Hz⁻¹, which is assumed to be a distinct constant for each occultation event. Thus, the SNR to be assumed in the simulations is related to the measured SNR as follows:

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}{\mathrm{SNR}}_{\mathrm{W},\mathrm{WOP}}& ={\displaystyle \frac{P_{\mathrm{signal}}}{P_{\mathrm{noise},,\mathrm{WOP}}}}={\displaystyle \frac{P_{\mathrm{signal}}}{P_{\mathrm{noise}}}}{\displaystyle \frac{B}{B_{\mathrm{WOP}}}}\\ & ={\mathrm{SNR}}_{\mathrm{W}}{\displaystyle \frac{B}{B_{\mathrm{WOP}}}}={\mathrm{SNR}}_{\mathrm{W}}{\displaystyle \frac{f_{\mathrm{s}}}{f_{\mathrm{s},\mathrm{WOP}}}}.\end{array}\end{array}$$

Then, the final formula for the noise amplitude to be added to the WOP signal is

$$\begin{array}{}\text{(A16)}& \begin{array}{rl}{\mathit{\sigma }}_{\mathrm{0}}{}^{\prime }& \approx \sqrt{{\displaystyle \frac{〈\widehat{u}{\widehat{u}}^{\ast }〉}{{\mathrm{SNR}}_{\mathrm{W},\mathrm{WOP}}}}}=\sqrt{{\displaystyle \frac{〈\widehat{u}{\widehat{u}}^{\ast }〉}{{\mathrm{SNR}}_{\mathrm{W}}}}{\displaystyle \frac{B_{\mathrm{WOP}}}{B}}}\\ & =\sqrt{{\displaystyle \frac{〈\widehat{u}{\widehat{u}}^{\ast }〉}{{\mathrm{SNR}}_{\mathrm{W}}}}{\displaystyle \frac{f_{\mathrm{s},\mathrm{WOP}}}{f_{\mathrm{s}}}}}.\end{array}\end{array}$$

Conventionally, the SNR is described in terms of voltage ratio in the RO community. In this way,

$$\begin{array}{}\text{(A17)}& {\mathrm{SNR}}_{\mathrm{W}}\left[\mathrm{W}{\mathrm{W}}^{-\mathrm{1}}\right]={\mathrm{SNR}}_V^{\mathrm{2}}\left[\mathrm{V}{\mathrm{V}}^{-\mathrm{1}}\right].\end{array}$$

Finally,

$$\begin{array}{}\text{(A18)}& {\mathit{\sigma }}_{\mathrm{0}}{}^{\prime }\approx \sqrt{{\displaystyle \frac{〈\widehat{u}{\widehat{u}}^{\ast }〉}{{\mathrm{SNR}}_V^{\mathrm{2}}}}{\displaystyle \frac{f_{\mathrm{s},\mathrm{WOP}}}{f_{\mathrm{s}}}}},\end{array}$$

which completes the derivation for the noise signal strength to be added to WOP signals.

The instrument noise added to WOP signals assumed the SNR of a MetOp-A occultation event to be the reference in Eq. (A18). The measurement is part of an end-of-life experimental campaign performed by EUMETSAT (European Organization for the Exploitation of Meteorological Satellites), where the vertical coverage of the GRAS (Global Navigation Satellite System Receiver for Atmospheric Sounding) instrument was extended temporarily up to 600 km SLTA (originally ∼ 80 km SLTA). Figure A1 shows the L1 C/A SNR of the occultation event scaled to f_s=1 Hz, which was not affected by ionospheric disturbances (S₄≤0.2), and the WOP amplitude with added noise on the last PS.

In our WOP simulations, the GNSS signal is propagated up to the rightmost PS. In order to define f_s,WOP in this particular scenario, the scanning velocity was approximated to v_s=3.2 km s⁻¹. Given the number of points per PS (2¹⁸) and the screen height (1000 km), the WOP sampling frequency in Eq. (A18) is $f_{\mathrm{s},\mathrm{WOP}}=\mathrm{839}$ Hz.

The S₄ index presented throughout the evaluations includes the added instrument noise. Thus, as in Syndergaard (2006),

$$\begin{array}{}\text{(A19)}& S_{\mathrm{4}}={\displaystyle \frac{\sqrt{〈{\left(I-〈\stackrel{\mathrm{‾}}{I}〉\right)}^{\mathrm{2}}〉}}{\stackrel{\mathrm{‾}}{I}}},\end{array}$$

where the signal intensity $I\propto \left|{\widehat{u}}_{\mathrm{total}}{\widehat{u}}_{\mathrm{total}}^{\ast }\right|$, $\stackrel{\mathrm{‾}}{I}$ stands for the filtered intensity, and 〈 〉 correspond to 10 s average.