the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Detection and localization of Flayer ionospheric irregularities with the backpropagation method along the radio occultation ray path
Vinícius LudwigBarbosa
Joel Rasch
Thomas Sievert
Anders Carlström
Mats I. Pettersson
Viet Thuy Vu
Jacob Christensen
The backpropagation (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 multipleirregularity 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 multiplebubble 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.
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 lowEarth 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 lowEarth 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 bendingangle inversion in the lower atmosphere (Gorbunov and Gurvich, 1998a, b; Dahl Mortensen, 1998). The regions with sharp gradients in refractivity, i.e. nonhomogeneities, 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 backpropagation (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 backpropagation 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 midlatitudes 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 postsunset 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 Fregion. 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 upflowing 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, solarcycle, and geomagneticactivity 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 Flayers 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 rescale 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 Flayer 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 Fregion 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.
Assuming the scenario of a GNSSRO 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):
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
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.
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.
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 backpropagation method, is written as in Gorbunov et al. (1996), Gorbunov and Gurvich (1998b), Dahl Mortensen (1998), and Sokolovskiy et al. (2002):
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
Thus, the final expression for Eq. (3) is
where $\mathrm{cos}\mathit{\xi}=\widehat{N}\cdot \widehat{\mathit{r}}$. Figure 3 shows the geometry of the backpropagation scenario, the relation between the angles, and the normal vector direction that changes along the LEO orbit.
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:

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

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 rightmost PS (via diffraction integral) by applying direct and inverse Fourier transforms recursively, viz.
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):
where $\mathit{u}{}^{\prime}=\mathit{u}\stackrel{\mathrm{\u203e}}{u}$ corresponds to detrended BP amplitude achieved by a threepass secondorder 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 secondaryfield 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.
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
where f is the carrier frequency, ρ is the electron density (el m^{−3}), 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^{−2} – 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 Fregion irregularity: plasma bubbles
Under low ionospheric activity, EDPs tend to resemble a slow function (Culverwell and Healy, 2015). Under highactivity 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 smallscale 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
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 \phantom{\rule{0.125em}{0ex}}\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 highorder terms composing the bias after the standard ionospheric correction (Vorob'ev and Krasil'nikova, 1994). The highorder bias, critical in meteorological and climate applications, is handled by either kappa or bilocal 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 backpropagation (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 veryhighfrequency (VHF) groundbased 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 2D random realization of Gaussian variables filtered by the spectral density function (SDF):
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 outerscale wave number, Γ is the Euler's gamma function, and ν denotes the spectral slope. The filtered variables,
are modulated to the electron density model,
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,
in which the function maximum and the bell width are set by
where x_{0} denotes the bubble placement along the x axis, hmF2 is the Fregion 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 LudwigBarbosa 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 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.
The simulated total field disturbed by the plasma bubble during the forward propagation, and sampled at the rightmost 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_{0} 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 singlebubble 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_{0} 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.
In the simulations, the filtered random field was modulated with an EDP modelled by Chapman's function (Culverwell and Healy, 2015) considering the Fregion peak ($nmF\mathrm{2}=\mathrm{8.81}\times {\mathrm{10}}^{\mathrm{11}}$ el m^{−3}), height (hmF2=288.5 km), and scale height (H=31 km) according to the EDP described in Carrano et al. (2011).
The forwardpropagation 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 (MetOpA) occultation event at a low latitude as a reference to the signaltonoise (SNR) level (see Fig. A1 in the Appendix). This measurement extends up to 600 km straightline 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 MetOpA endoflife 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 Fregion peak was around 600 V V^{−1}.
In forwardpropagation simulations, the closest phase screen to x_{0}, 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.
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}\phantom{\rule{0.125em}{0ex}}\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}}\phantom{\rule{0.125em}{0ex}}{\mathrm{el}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{\mathrm{3}}$, which results in weak scattering (S_{4}=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},\mathrm{min}}=\mathrm{1.7}$ km.
Figure 7 shows the result considering the bubble placement on the outbound ray path.
In the singlebubble scenario, the location estimate has good accuracy regardless of the placement in the inbound or outbound sector. Therefore, the location estimate in a singlebubble scenario is limited by the precision considered in the BP method, herein 5 km. The minor difference in the scintillation index (S_{4}) 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.
The curve corresponding to ${\mathit{\sigma}}_{\mathrm{\Delta}\mathit{\rho}/\mathit{\rho}}\le \mathrm{2}\phantom{\rule{0.125em}{0ex}}\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}≈0.0456 after Figs. 8b and A1b). Thus, this indicates that the estimations are unreliable when σ_{u}≤σ_{0}. For rms fluctuation levels ${\mathit{\sigma}}_{\mathrm{\Delta}\mathit{\rho}/\mathit{\rho}}>\mathrm{3.0}\phantom{\rule{0.125em}{0ex}}\mathit{\%}$ ($\pm \mathrm{2.64}\times {\mathrm{10}}^{\mathrm{10}}$ el m^{−3}), 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}<\phantom{\rule{0.125em}{0ex}}\mathrm{146.7}$ km). Regardless, the disturbance level ${\mathit{\sigma}}_{\mathrm{\Delta}\mathit{\rho}/\mathit{\rho}}=\mathrm{3.0}\phantom{\rule{0.125em}{0ex}}\mathit{\%}$ assumed in the forward propagation corresponds to a weal disturbance in the LEO's orbit (S_{4}<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}\phantom{\rule{0.125em}{0ex}}\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).
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}\phantom{\rule{0.125em}{0ex}}\mathit{\%}$). A region with extension L_{H}≤20 km creates low scintillation in the GNSS signal (S_{4}<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.
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).
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}\phantom{\rule{0.125em}{0ex}}\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.
A comparison between Figs. 10c, d (widebubble 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 singleregion cases where the predominant constraint to detection is the noise level (σ_{u}≈σ_{0}), these results indicate that the separation between the regions has a major influence on the detection and location task of multiple patches.
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}\phantom{\rule{0.125em}{0ex}}\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}\phantom{\rule{0.125em}{0ex}}\mathit{\%}$ and ${\mathit{\sigma}}_{\mathrm{\Delta}\mathit{\rho}/\mathit{\rho}}=\mathrm{12}\phantom{\rule{0.125em}{0ex}}\mathit{\%}$.
The standard deviation curves in scenarios including a bubble with ${\mathit{\sigma}}_{\mathrm{\Delta}\mathit{\rho}/\mathit{\rho}}=\mathrm{6}\phantom{\rule{0.125em}{0ex}}\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.
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 sporadicE scintillation (Zeng and Sokolovskiy, 2010; Arras and Wickert, 2018; Yu et al., 2020; Carmona et al., 2022). In Fig. 15b, a ushaped 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 ushaped 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).
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 singlebubble 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}\phantom{\rule{0.125em}{0ex}}\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_{4}>0.2. However, the accuracy of the centre estimate decreases with increasing width.
In multiplebubble 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 forwardpropagation 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}<σ_{0}). 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^{−3}). 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 MetOpA 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 groundbased 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.
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.
The total signal will also include noise,
We used the measured SNR from a representative MetOpA occultation event to estimate the appropriate noise level added to the WOP amplitude,
The noise in occultation measurements has several sources: thermal noise in the receiver, clock noise, and cochannel noise. For this task, we assumed a normal distribution to model the white noise – i.e. $X,Y\phantom{\rule{0.125em}{0ex}}\sim \mathcal{N}(\mathit{\mu},{\mathit{\sigma}}^{\mathrm{2}})$, where μ is the mean value, and σ^{2} is the variance. Then, the noise in the ith sample is
where $X,Y\phantom{\rule{0.125em}{0ex}}\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,
Likewise, the averaged signal power becomes
The SNR in terms of the signal and noise power, with units [W W^{−1}], is given by
Hence,
and
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
in which f_{s} is the sample rate in Hz. The noise power is given by
where N_{0} is the noise power density in W Hz^{−1}, 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:
Then, the final formula for the noise amplitude to be added to the WOP signal is
Conventionally, the SNR is described in terms of voltage ratio in the RO community. In this way,
Finally,
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 MetOpA occultation event to be the reference in Eq. (A18). The measurement is part of an endoflife 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_{4}≤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^{−1}. Given the number of points per PS (2^{18}) and the screen height (1000 km), the WOP sampling frequency in Eq. (A18) is ${f}_{\mathrm{s},\mathrm{WOP}}=\mathrm{839}$ Hz.
The S_{4} index presented throughout the evaluations includes the added instrument noise. Thus, as in Syndergaard (2006),
where the signal intensity $I\propto \left{\widehat{u}}_{\mathrm{total}}\phantom{\rule{0.125em}{0ex}}{\widehat{u}}_{\mathrm{total}}^{\ast}\right$, $\stackrel{\mathrm{\u203e}}{I}$ stands for the filtered intensity, and 〈 〉 correspond to 10 s average.
Codes and measurements are provided by request to the authors.
VLB, JR, and TS designed the study cases, and VLB performed the simulations and processing. VLB prepared the paper and its revised versions with contributions from all the coauthors. JR prepared the Appendix.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to thank Riccardo Notarpietro (EUMETSAT) for sharing the GRAS ionospheric extension experiment on EPS MetOpA data, which was used as a reference in this work.
This research has been supported by the Swedish National Space Agency (NRFP4 program).
This paper was edited by Peter Alexander and reviewed by two anonymous referees.
Aarons, J.: Global morphology of ionospheric scintillations, Proceedings of the IEEE, 70, 360–378, https://doi.org/10.1109/PROC.1982.12314, 1982. a, b
Abdu, M. A., Nogueira, P. A. B., Santos, A. M., de Souza, J. R., Batista, I. S., and Sobral, J. H. A.: Impact of disturbance electric fields in the evening on prereversal vertical drift and spread F developments in the equatorial ionosphere, Ann. Geophys., 36, 609–620, https://doi.org/10.5194/angeo366092018, 2018. a
Arras, C. and Wickert, J.: Estimation of ionospheric sporadic E intensities from GPS radio occultation measurements, J. Atmos. Sol.Terr. Phys., 171, 60–63, https://doi.org/10.1016/j.jastp.2017.08.006, 2018. a, b, c
Béniguel, Y., Romano, V., Alfonsi, L., Aquino, M., Bourdillon, A., Cannon, P., Franceschi, G. D., Dubey, S., Forte, B., Gherm, V., Jakowski, N., Materassi, M., Noack, T., Pozoga, M., Rogers, N., Spalla, P., Strangeways, H. J., Warrington, E. M., Wernik, A. W., Wilken, V., and Zernov, N.: Ionospheric scintillation monitoring and modelling, Ann. Geophys., 52, 391–416, https://doi.org/10.4401/ag4595, 2009. a
Bevis, M., Businger, S., Herring, T. A., Rocken, C., Anthes, R. A., and Ware, R. H.: GPS meteorology: Remote sensing of atmospheric water vapor using the global positioning system, J. Geophys. Res., 97, 15787, https://doi.org/10.1029/92JD01517, 1992. a
Carmona, R. A., Nava, O. A., Dao, E. V., and Emmons, D. J.: A Comparison of SporadicE Occurrence Rates Using GPS Radio Occultation and Ionosonde Measurements, Remote Sens., 14, 1–15, https://doi.org/10.3390/rs14030581, 2022. a, b
Carrano, C. S., Groves, K. M., Caton, R. G., Rino, C. L., and Straus, P. R.: Multiple phase screen modeling of ionospheric scintillation along radio occultation raypaths, Radio Sci., 46, 1–14, https://doi.org/10.1029/2010RS004591, 2011. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o
Carrano, C. S., Groves, K. M., Mcneil, W. J., and Doherty, P. H.: Scintillation Characteristics across the GPS Frequency Band, in: 25th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS), 17–21 September 2012, Nashville, TN, Institute of Navigation, 1972–1989, ISBN: 9781622769803, 2012. a
Carrano, C. S., Groves, K. M., Delay, S. H., and Doherty, P. H.: An inverse diffraction technique for scaling measurements of ionospheric scintillations on the GPS L1, L2, and L5 carriers to other frequencies, in: Institute of Navigation International Technical Meeting (ITM), 27–29 January, San Diego, CA, Institute of Navigation, 709–719, ISSN: 23303646, 2014. a
Cherniak, I. and Zakharenkova, I.: Highlatitude ionospheric irregularities: differences between ground and spacebased GPS measurements during the 2015 St. Patrick’s Day storm, Earth Planet. Space, 68, 136, https://doi.org/10.1186/s4062301605061, 2016. a
Cherniak, I., Zakharenkova, I., and Sokolovskiy, S.: Multi‐Instrumental Observation of Storm‐Induced Ionospheric Plasma Bubbles at Equatorial and Middle Latitudes, J. Geophys. Res.Space, 124, 1491–1508, https://doi.org/10.1029/2018JA026309, 2019. a, b, c, d, e, f, g
Culverwell, I. D. and Healy, S. B.: Simulation of L1 and L2 bending angles with a model ionosphere, Tech. Rep. 17, Danish Meterological Institute, Copenhagen, http://www.romsaf.org/generaldocuments/rsr/rsr_17.pdf (last access: 30 March 2023), 2015. a, b, c, d
Dahl Mortensen, M.: The BackPropagation Method for Inversion of Radio Occultation Data, Tech. Rep. 14, Danish Meteorological Institute, Copenhagen, Denmark, https://www.dmi.dk/fileadmin/Rapporter/SR/sr9814.pdf (last access: 30 March 2023), 1998. a, b
Fjeldbo, G., Kliore, A. J., and Eshleman, V. R.: The Neutral Atmosphere of Venus as Studied with the Mariner V Radio Occultation Experiments, Astron. J., 76, 123–140, https://doi.org/10.1086/111096, 1971. a
Gorbunov, M. and Gurvich, A. S.: Microlab1 experiment: Multipath effects in the lower troposphere, J. Geophys. Res.Atmos., 103, 13819–13826, https://doi.org/10.1029/98JD00806, 1998a. a, b
Gorbunov, M. and Gurvich, A. S.: Algorithms of inversion of microlab1 satellite data including effects of multipath propagation, Int. J. Remote Sens., 19, 2283–2300, https://doi.org/10.1080/014311698214721, 1998b. a, b
Gorbunov, M. E. and Lauritsen, K. B.: Linearized Zverev Transform and its application for modeling radio occultations, Radio Sci., 42, RS3023, https://doi.org/10.1029/2006RS003590, 2007. a, b
Gorbunov, M. E., Gurvich, A. S., and Bengtsson, L.: Advanced Algorithms of Inversion of GPS/MET Satellite Data and Their Application to Reconstruction of Temperature and Humidity, Tech. Rep. 211, MaxPlanck Institute for Meteorology, Hamburg, Germany, https://pure.mpg.de/pubman/faces/ViewItemOverviewPage.jsp?itemId=item_3033351 (last access: 30 March 2023), 1996. a, b
Gorbunov, M. E., Gurvich, A. S., and Shmakov, A. V.: Backpropagation and radioholographic methods for investigation of sporadic ionospheric Elayers from Microlab1 data, Int. J. Remote Sens., 23, 675–685, https://doi.org/10.1080/01431160010030091, 2002. a, b, c, d, e
Healy, S. B. and Culverwell, I. D.: A modification to the standard ionospheric correction method used in GPS radio occultation, Atmos. Meas. Tech., 8, 3385–3393, https://doi.org/10.5194/amt833852015, 2015. a
Jiao, Y. and Morton, Y. T.: Comparison of the effect of high‐latitude and equatorial ionospheric scintillation on GPS signals during the maximum of solar cycle 24, Radio Sci., 50, 886–903, https://doi.org/10.1002/2015RS005719, 2015. a, b
Kelley, M. C., Larsen, M. F., LaHoz, C., and McClure, J. P.: Gravity wave initiation of equatorial spread F: A case study, J. Geophys. Res., 86, 9087, https://doi.org/10.1029/JA086iA11p09087, 1981. a
Kelly, M. A., Comberiate, J. M., Miller, E. S., and Paxton, L. J.: Progress toward forecasting of space weather effects on UHF SATCOM after Operation Anaconda, Space Weather, 12, 601–611, https://doi.org/10.1002/2014SW001081, 2014. a
Kepkar, A., Arras, C., Wickert, J., Schuh, H., Alizadeh, M., and Tsai, L.C.: Occurrence climatology of equatorial plasma bubbles derived using FormoSat3/COSMIC GPS radio occultation data, Ann. Geophys., 38, 611–623, https://doi.org/10.5194/angeo386112020, 2020. a, b
Knepp, D.: Multiple phasescreen calculation of the temporal behavior of stochastic waves, Proceedings of the IEEE, 71, 722–737, https://doi.org/10.1109/PROC.1983.12660, 1983. a, b
Kursinski, E. R., Hajj, G. A., Schofield, J. T., Linfield, R. P., and Hardy, K. R.: Observing Earth's atmosphere with radio occultation measurements using the Global Positioning System, J. Geophys. Res.Atmos., 102, 23429–23465, https://doi.org/10.1029/97JD01569, 1997. a, b
Liu, C., Kirchengast, G., Syndergaard, S., Schwaerz, M., Danzer, J., and Sun, Y.: New higher‐order correction of GNSSRO bending angles accounting for ionospheric asymmetry: Evaluation of performance and added value, Remote Sens., 12, 1–23, https://doi.org/10.3390/rs12213637, 2020. a
LudwigBarbosa, V., Rasch, J., Carlström, A., Pettersson, M. I., and Vu, V. T.: GNSS Radio Occultation Simulation Using Multiple Phase Screen Orbit Sampling, IEEE Geosci. Remote Sens. Lett., 17, 1323–1327, https://doi.org/10.1109/LGRS.2019.2944537, 2020. a
Ludwig‐Barbosa, V., Sievert, T., Rasch, J., Carlström, A., Pettersson, M. I., and Thuy Vu, V.: Evaluation of Ionospheric Scintillation in GNSS Radio Occultation Measurements and Simulations, Radio Sci., 55, e2019RS006996, https://doi.org/10.1029/2019RS006996, 2020. a
Ma, G., Hocke, K., Li, J., Wan, Q., Lu, W., and Fu, W.: GNSS Ionosphere Sounding of Equatorial Plasma Bubbles, Atmosphere, 10, 1–11, https://doi.org/10.3390/atmos10110676, 2019. a
Resende, L. C. A., Arras, C., Batista, I. S., Denardini, C. M., Bertollotto, T. O., and Moro, J.: Study of sporadic E layers based on GPS radio occultation measurements and digisonde data over the Brazilian region, Ann. Geophys., 36, 587–593, https://doi.org/10.5194/angeo365872018, 2018. a
Sokolovskiy, S., Schreiner, W., Rocken, C., and Hunt, D.: Detection of highaltitude ionospheric irregularities with GPS/MET, Geophys. Res. Lett., 29, 1–4, https://doi.org/10.1029/2001GL013398, 2002. a, b, c, d, e, f, g
Sokolovskiy, S., Schreiner, W. S., Zeng, Z., Hunt, D. C., Kuo, Y.H., Meehan, T. K., Stecheson, T. W., Mannucci, A. J., and Ao, C. O.: Use of the L2C signal for inversions of GPS radio occultation data in the neutral atmosphere, GPS Solutions, 18, 405–416, https://doi.org/10.1007/s102910130340x, 2014. a, b, c
Sommerfeld, A.: Optics: Lectures on Theoretical Physics, Vol. IV, Academic Press Inc., New York, USA, ISBN 0126546762, 1964. a, b, c
Stolle, C., Lühr, H., Rother, M., and Balasis, G.: Magnetic signatures of equatorial spread F as observed by the CHAMP satellite, J. Geophys. Res., 111, A02304, https://doi.org/10.1029/2005JA011184, 2006. a, b
Stolle, C., Lühr, H., and Fejer, B. G.: Relation between the occurrence rate of ESF and the equatorial vertical plasma drift velocity at sunset derived from global observations, Ann. Geophys., 26, 3979–3988, https://doi.org/10.5194/angeo2639792008, 2008. a
Syndergaard, S.: COSMIC S4 Data, https://cdaacwww.cosmic.ucar.edu/cdaac/doc/documents/s4_description.pdf (last access: 30 March 2023), 2006. a
Vorob'ev, V. V. and Krasil'nikova, T. G.: Estimation of the accuracy of the atmospheric refractive index recovery from doppler shift measurements at frequencies used in the NAVSTAR system, USSR Phys. Atmos. Ocean, Engl. Transl., 29, 602–609, 1994. a
Wickert, J., Pavelyev, A. G., Liou, Y. A., Schmidt, T., Reigber, C., Igarashi, K., Pavelyev, A. A., and Matyugov, S. S.: Amplitude variations in GPS signals as a possible indicator of ionospheric structures, Geophys. Res. Lett., 31, L24801, https://doi.org/10.1029/2004GL020607, 2004. a, b
Xiong, C., Stolle, C., Lühr, H., Park, J., Fejer, B. G., and Kervalishvili, G. N.: Scale analysis of equatorial plasma irregularities derived from Swarm constellation, Earth Planet. Space, 68, 121, https://doi.org/10.1186/s4062301605025, 2016. a, b
Yeh, K. and Liu, C.H.: Radio wave scintillations in the ionosphere, Proceedings of the IEEE, 70, 324–360, https://doi.org/10.1109/PROC.1982.12313, 1982. a, b, c
Yu, B., Xue, X., Yue, X., Yang, C., Yu, C., Dou, X., Ning, B., and Hu, L.: The global climatology of the intensity of the ionospheric sporadic E layer, Atmos. Chem. Phys., 19, 4139–4151, https://doi.org/10.5194/acp1941392019, 2019. a
Yu, B., Scott, C. J., Xue, X., Yue, X., and Dou, X.: Derivation of global ionospheric Sporadic E critical frequency f0Es data from the amplitude variations in GPS/GNSS radio occultations, Roy. Soc. Open Sci., 7, 200320, https://doi.org/10.1098/rsos.200320, 2020. a, b
Zeng, Z. and Sokolovskiy, S.: Effect of sporadic E clouds on GPS radio occultation signals, Geophys. Res. Lett., 37, 1–5, https://doi.org/10.1029/2010GL044561, 2010. a, b, c
Zeng, Z., Sokolovskiy, S., Schreiner, W. S., and Hunt, D.: Representation of Vertical Atmospheric Structures by Radio Occultation Observations in the Upper Troposphere and Lower Stratosphere: Comparison to HighResolution Radiosonde Profiles, J. Atmos. Ocean. Technol., 36, 655–670, https://doi.org/10.1175/JTECHD180105.1, 2019. a