the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Analysis of ionospheric structure influences on residual ionospheric errors in GNSS radio occultation bending angles based on ray tracing simulations
Congliang Liu
Gottfried Kirchengast
Yueqiang Sun
Kefei Zhang
Robert Norman
Marc Schwaerz
Weihua Bai
Qifei Du
Ying Li
The Global Navigation Satellite System (GNSS) radio occultation (RO) technique is widely used to observe the atmosphere for applications such as numerical weather prediction and global climate monitoring. The ionosphere is a major error source to RO at upper stratospheric altitudes, and a linear dualfrequency bending angle correction is commonly used to remove the firstorder ionospheric effect. However, the higherorder residual ionospheric error (RIE) can still be significant, so it needs to be further mitigated for highaccuracy applications, especially from 35 km altitude upward, where the RIE is most relevant compared to the decreasing magnitude of the atmospheric bending angle. In a previous study we quantified RIEs using an ensemble of about 700 quasirealistic endtoend simulated RO events, finding typical RIEs at the 0.1 to 0.5 µrad noise level, but were left with 26 exceptional events with anomalous RIEs at the 1 to 10 µrad level that remained unexplained. In this study, we focused on investigating the causes of the high RIE of these exceptional events, employing detailed alongraypath analyses of atmospheric and ionospheric refractivities, impact parameter changes, and bending angles and RIEs under asymmetric and symmetric ionospheric structures. We found that the main causes of the high RIEs are a combination of physicsbased effects – where asymmetric ionospheric conditions play the primary role, more than the ionization level driven by solar activity – and technical ray tracer effects due to occasions of imperfect smoothness in ionospheric refractivity model derivatives. We also found that alongray impact parameter variations of more than 10 to 20 m are possible due to ionospheric asymmetries and, depending on prevailing horizontal refractivity gradients, are positive or negative relative to the initial impact parameter at the GNSS transmitter. Furthermore, mesospheric RIEs are found generally higher than upperstratospheric ones, likely due to being closer in tangent point heights to the ionospheric E layer peaking near 105 km, which increases RIE vulnerability. In the future we will further improve the alongray modeling system to fully isolate technical from physicsbased effects and to use it beyond this work for additional GNSS RO signal propagation studies.
 Article
(10617 KB)  Companion paper
 BibTeX
 EndNote
Global Navigation Satellite System (GNSS) radio occultation (RO; Melbourne et al., 1994; Kursinski et al., 1997; Hajj et al., 2002) is a relatively new atmospheric sounding technique. It can deliver data traceable to the international standard of time (the SI second) and has a demonstrated capacity for monitoring decadalscale climate change in the free atmosphere (Steiner et al., 2009, 2011, 2013; Anthes, 2011; Foelsche et al., 2011; Lackner et al., 2011; Ho et al., 2012; Angerer et al., 2017). This capacity rests on RO's unique combination of characteristics such as high vertical resolution, high accuracy, longterm stability, and global coverage (Kursinski et al., 1997; ScherllinPirscher et al., 2011; Anthes, 2011; Steiner et al., 2011). Figure 1 illustrates the GNSS RO geometry that constitutes the basis of the RO technique. The focus is to schematically show essential aspects relevant to this study on alongray ionospheric influences on RO bending angles, which deepens insight on top of our recent Liu et al. (2015) study.
Ionospheric error is significant in GNSS RO observations (e.g., Kursinski et al., 1997; Mannucci et al., 2011; Liu et al., 2013), and a dualfrequency linear combination of RO bending angles is usually implemented to correct for the firstorder ionospheric effect (Vorob'ev and Krasil'nikova, 1994; Ladreiter and Kirchengast, 1996). However, the higherorder residual ionospheric error (RIE) after this correction is still not negligible for highaccuracy applications such as RObased climate change monitoring (Steiner et al., 2011, 2013). This applies especially above about 35 km altitude, where the RIE becomes increasingly relevant compared to the exponentially decreasing magnitude of the neutral atmospheric bending angle (Syndergaard, 2000; Mannucci et al., 2011; Danzer et al., 2013, 2015; Liu et al., 2013, 2015; Healy and Culverwell, 2015).
Moreover, the RIE can propagate downwards into the lowerstratospheric retrievals of refractivity and temperature through the Abel integral and the hydrostatic integral (Kursinski et al., 1997; Gobiet and Kirchengast, 2004; Steiner and Kirchengast, 2005; Gobiet et al., 2007). It is therefore essential to better understand and further mitigate the RIE in order to enable benchmarkquality stratospheric RO retrievals.
A wide array of studies related to a better understanding of higherorder ionospheric errors in GNSS RO data have been conducted already by a range of scientists (Bassiri and Hajj, 1993; Vorob'ev and Krasil'nikova, 1994; Ladreiter and Kirchengast, 1996; Syndergaard, 2000; Gorbunov, 2002; Hoque and Jakowski, 2010, 2011; Mannucci et al., 2011; Danzer et al., 2013, 2015; Healy and Culverwell, 2015; Coleman and Forte, 2017). A few of these also suggested ways of correcting higherorder RIEs in RO bending angles (Syndergaard, 2000; Danzer et al., 2013; Healy and Culverwell, 2015), which may be applied on top of the standard dualfrequency correction introduced by Vorob'ev and Krasil'nikova (1994).
The convenient formulation introduced by Healy and Culverwell (2015), which adds a fairly simple higherorder squaredbendingangle difference term to the standard correction, is meanwhile applied in operational processing of the data from the European MetOp (Meteorological Operational Satellites) RO mission (Luntama et al., 2008; Christian Marquardt, EUMETSAT Darmstadt, personal communications, 2017). Recently, Angling et al. (2018) further improved the empirical modeling of the “kappa coefficient” in this formulation, by accounting for solar zenith angle, solar flux (F10.7 index), and altitude dependencies.
In our work over the recent years we have assessed the variation of bending angle RIEs (biases and standard deviations) with solar activity, with latitudinal region, and with or without the assumption of ionospheric spherical symmetry and of coexisting RO observing system errors, using endtoend simulations for single RO events (Liu et al., 2013) and a fullday ensemble of RO events (Liu et al., 2015). As shown in these explanatory simulation studies, in overall agreement with the empirical study of Danzer et al. (2013), the RIE biases have a clear negative tendency and a bias magnitude increasing with solar activity, as well as being affected by deviations from ionospheric spherical symmetry (Mannucci et al., 2010) where increasing asymmetries also tend to increase the biases.
What remained unexplored in our Liu et al. (2015) study and had also not yet been explored elsewhere – but is critical to be understood for further improvement of the existing RIE corrections that apply spherical symmetry (Syndergaard, 2000; Healy and Culvervell, 2015; Angling et al., 2018) – is the influences of the threedimensional and asymmetric ionospheric structures along the GNSStoLEO (low Earth orbit) signal paths on the RIE, in particular the conditions that may lead to anomalously high RIEs.
A first step in this direction, though not focusing on bending angle RIEs, was the study by Mannucci et al. (2011), which found that under ionospheric storm conditions anomalous effects can be significant. Recently also Coleman and Forte (2017) reported RIE investigations for asymmetry conditions, including on the effect of traveling ionospheric disturbances upon the RIE. Another step was the somewhat puzzling side result in our Liu et al. (2015) study that the endtoend simulations of an ensemble of about 700 RO events produced about two dozen RIE outlier profiles. The basis was 3D ray tracing simulations, where the ionospheric model NeUoG (University of Graz electron density model; Leitinger and Kirchengast, 1997) was used as a quasirealistic model for largescale 3D ionospheric structures, together with the atmospheric model MSIS90 (Mass Spectrometer and Incoherent Scatter neutral atmosphere model 1990; Hedin, 1991) for simple but representative neutral atmosphere reference conditions. More precisely, the RIE standard deviation of 26 profiles from the simulations exceeded a threshold value of 7 µrad within the upper stratosphere and mesosphere. These were therefore rejected from the ensemble statistics results reported by Liu et al. (2015).
In this study we now place focus on these 26 exceptional profiles and, by way of detailed alongray analyses of ray tracing simulations, aim to shed light on the causes of anomalously high RIEs, with the additional goal of deepening quantitative insight into how RIEs accumulate during signal propagation, along with accumulation of the total (atmospheric) bending angles that are the desired RO observables. In Sect. 2, the exceptional RO events and the simulation setup for exploring their bending angle RIEs are introduced. Section 3 provides the results, which we mainly discuss through detailed inspection of example events. A summary and conclusions are finally given in Sect. 4.
2.1 Exceptional RO events
The ensemble of RO events used by Liu et al. (2015) was simulated for 15 July 2008, adopting the European MetOp RO mission as an example low Earth orbiter (Edwards and Pawlak, 2000), specifically thinking of MetOpA, which was launched as the first of the MetOp series in late 2006 (Luntama et al., 2008). Each MetOp satellite is a sunsynchronous LEO satellite at about 820 km with the Global Positioning System (GPS) Receiver for Atmospheric Sounding (GRAS) on board (Loiselet et al., 2000), which acquires about 700 RO events per day (Luntama et al., 2008).
Using, as summarized above, simple spherically symmetric neutral atmospheric modeling (by MSIS90) combined with 3D ionospheric modeling (by NeUoG), we simulated in that study the ensemble of daily RO events for 14 different endtoend simulation cases. These included withoutionosphere (wi) cases as well as spherical symmetry (ss) and nonsphericalsymmetry (ns) ionospheric cases for low, medium, and high solar activity levels, under the assumption of either perfect observing system (op) with no errors or realistic observing system (or) with MetOptype errors; for details see Liu et al. (2015), Table 2 and Sect. 2.3 therein. The total number of the simulated RO events found for the day was 723, of which 26 exceptionally noisy ones were classified as outliers (estimated bending angle RIE exceeding 7 µrad somewhere within 30 to 80 km). These 26 events are investigated closer in this study.
Figure 2a shows the global distribution of mean tangent point (TP) locations of all 723 events (as small triangles) and highlights the locations of the 26 exceptional events (as red triangles). The majority of the latter (18 of the 26) appear to cluster over the European–Asian and Indian Ocean regions (EAC and IOC, highlighted as boxes); the remaining eight events are distributed more diversely in other extratropical locations, mainly in the Northern Hemisphere. Figure 2b and c depict the RIE bias and standard deviation, defined in the same way as by Liu et al. (2013), which are estimated for the upper stratosphere and mesosphere (30–80 km) for the 26 events, for the nonsphericalsymmetry (“opns”) and spherical symmetry (“opss”) ionospheric conditions, respectively. Intercomparing Fig. 2b and c shows that the main driver of anomalously high RIEs is asymmetric ionospheric conditions and possibly residual error effects from ray tracing through the 3D ionosphere, since only few events (6 of the 26) exhibit large RIE standard deviations (exceeding 1 µrad) even in the case of symmetric ionospheric conditions.
Related to the clusters, one can see that, in the opss case, almost all noisy exceptional events occurred in the IOC, while in the opns case the noisiest ones occurred in both the EAC and IOC. Related to solar activity, one can see that in both the opss and opns cases higher ionization (F10.7) levels generally lead to increased RIEs, compared to lowest ionization (F10.7 = 70), but the picture is ambiguous, and often medium solar activity also leads to higher RIEs than high solar activity.
These overall characteristics revealed by Fig. 2 point, in particular, to two facts that shall guide our detailed investigation for better understanding of anomalous RIEs: (1) asymmetric ionospheric conditions play a key role, more than ionization levels and possible geographic location dependencies (e.g., via solar or geomagnetic influences), and so inspection of the alongray signal dynamics is essential; (2) the several exceptions from the overall characteristics, and some geographic clustering that has no obvious physicsbased explanation, indicate that there is no single clear cause for the anomalous RIEs and that some perturbations also come in from the technical challenge of smooth ray tracing at millimetric excessphase accuracy through 3D ionospheric models like NeUoG.
We inspected the bending angle RIE profiles of the 26 events over the 20 to 80 km height range, including also their underlying excessphase RIE profiles, and chose three representative events that we will explore in detail below for improving RIE insight: an extremely noisy event (Occ.530 from the EAC) and a medium noisy event (Occ.20 from the IOC) from the 26 exceptional events, both used at medium solar activity, and a reference event from the 697 standard events, with lownoise RIE (Occ.25). Table 1 summarizes the main parameters for these three events, and Fig. 3 illustrates them in terms of excess phases, bending angles, and the associated RIEs.
Figure 3a shows the behavior of the excess phases of the three events. The L1 and L2 excess phases are around −11 to −15 and −18 to −25 m, respectively, typical for medium solar activity (Liu et al., 2013). After the standard ionospheric correction, the ionospherecorrected (Lc) excess phases are found near 0 m as should be the case. The excessphase RIE profiles (Fig. 3b) exhibit some spiky behavior for the two exceptional events, on top of comparatively lownoise RIEs otherwise. This points to unphysical values at the spiky impact height levels, given that the largescale 3D ionospheric structure of the NeUoG model should be physically unable to induce such sharp changes. It hence indicates that the ray tracing is technically challenged along the signal propagation paths pertaining to these levels by slight ionospheric model discontinuities, which render millimetric excessphase accuracy unattainable for these ray paths.
As Leitinger and Kirchengast (1997) describe, substantial empirical modeling effort went into strict smoothness of the NeUoG electron density field and its 3D derivatives that are key for highaccuracy ray tracing; nevertheless some slight discontinuities have likely remained in some rare locations of the modeling space spanned by altitude, latitude, longitude, (universal) time, month, and solar activity. It will therefore be important to separate such technical modeling effects from the physical effects on the propagating signals that cause high RIEs.
Figure 3c shows that, for all three events, the difference between α_{1} and α_{2} somewhat increases with increasing impact height, a feature already visible in the Liu et al. (2013) results. It is caused by the increasing ionospheric influence when tangent point heights gradually approach ionospheric E layer heights around 105 km from below. These overall differences between α_{1} and α_{2} amount to about 15 to 20 µrad near 80 km and are effectively eliminated by the standard ionospheric correction, bringing the α_{c} profile to near zero as should be the case. Nevertheless, substantial waveformlike perturbations remain on α_{c} for the two exceptional events Occ.530 and Occ.20, which even more clearly show up in the bending angle RIE profile (Fig. 3d).
Intercomparing Fig. 3d with b suggests that these waveformlike perturbations in the bending angle RIE are mainly induced by propagating the spiky excessphase perturbations through the bending angle retrieval, which involves filtering and a derivative operation from excess phase to Doppler shift (Schwarz et al., 2017). One main cause that has driven many of the exceptional events into the outlier range (i.e., into exceeding 7 µrad somewhere within 30 to 80 km) is thus evidently the technical effects from the ray tracing through the NeUoG ionosphere, which is not perfectly smooth everywhere in its electron density and hence refractivity field derivatives. It is thus important to more closely explore the alongray signal dynamics in order to understand how such technical effects may occur along ray paths and in particular in order to better understand the physical effects that drive high RIEs. Our related alongray analysis methodology is introduced next.
2.2 Investigation methodology
2.2.1 Ray tracing method
The ray tracing technique is commonly used for calculating propagation paths of an electromagnetic signal in a medium specified by a positiondependent refractive index field. It has become a significant tool for investigating signal propagation in RO technology. For example, Ladreiter and Kirchengast (1996), Syndergaard (2000), Gobiet and Kirchengast (2004), Steiner and Kirchengast (2005), Hoque and Jakowski (2010), Mannucci et al. (2011), Danzer et al. (2013, 2015), and Liu et al. (2013, 2015) have employed this method inter alia or with a main topical focus to investigate the ionospheric effects on GNSS RO signals. Danzer et al. (2015) noted that their analysis was somewhat limited by high noise of the simulated bending angle profiles at mid to high latitudes, which partly reflected the degrading impact of technical ray tracer effects that we also encounter and more explicitly address in this study.
We employ the 3D numerical ray tracing technique integrated in the Endtoend GNSS Occultation Performance Simulation and Processing System version 5.6 (EGOPS 5.6; Fritzer et al., 2013) in the same way as used by Liu et al. (2013, 2015) for simulating the GPStoLEO signal propagation through the atmosphere–ionosphere system; for a detailed description of the endtoend simulation setup the reader is therefore referred to these recent studies. Here we specifically refined and enhanced this setup in the 3D ray tracing part by adding the cocomputation and result extraction for a range of key variables along the propagation paths, instead of only providing the final observational variables of an RO event at the LEO receiver position.
2.2.2 Investigated variables
We implemented detailed alongray diagnostic capabilities into the 3D ray tracer of the EGOPS 5.6 software (Fritzer et al., 2013), which is an extensively proven highaccuracy ray tracer originally developed in the 1990s (Syndergaard, 1998, 1999). In particular, we computed the following key diagnostic variables for all individual numerical steps along the ray paths simulated for the GPS L1 and L2 frequencies as well as for a reference case without ionosphere (Lr), with each ray path starting at the GPS transmitter position and ending at the LEO receiver position: 3D position in the WGS84based Earthcentered, Earthfixed (ECEF) system, storing both the cartesian (X, Y, Z) and geodetic (height, latitude, longitude) coordinates; alongray distance relative to the TP, the latter evaluated as the ray's point of closest approach to the WGS84 ellipsoidal surface (parabolic vertex fit to the three alongray positions closest to the surface); atmospheric refractivity; L1 and L2 ionospheric refractivity; L1 and L2 impact parameter and impact parameter difference to the initial impact parameter at the GPS transmitter position (termed “delta impact parameter”, induced along the ray in the case of nonsphericalsymmetry conditions); accumulated L1, L2, and ionospherecorrected (Lc) bending angle (bending angle accrued from the GPS transmitter position to the alongray position); and RIE of the Lc bending angle, estimated relative to the Lr bending angle obtained from a simulation case without ionosphere (Liu et al., 2013).
These alongray variables are computed for all available ray paths from 80 to 20 km impact height, which are produced at 50 Hz sampling rate for any RO event investigated. This leads to a dense sampling by roughly 1500 ray paths in this altitude range (i.e., typical average scan velocities of RO events are near 2 km s^{−1} in this domain). Likewise the ray tracer provides fairly dense alongray stepping, employing an adaptive step size concept with finest steps at highest local refractive index changes (for details see, e.g., Syndergaard, 1999; Fritzer et al., 2013). Together these features enable inspecting the propagation characteristics of RO events through the atmosphere–ionosphere system at very high resolution in a convenient 2D alongray distance vs. impact height coordinate system that accurately represents the real 3Dwarped occultation event plane between the GPS and LEO orbit arcs.
We will inspect the results for the three representative RO events chosen (Occ.530, Occ.20, Occ.25; see Sect. 2.1 above) in this alongray distance vs. impact height coordinate system. Before turning to this, we briefly summarize here the equations for the alongray computation of those key variables that we will inspect closely. This aims to facilitate an appropriate understanding and interpretation of the results.
On the basis of Snell's law, when the Earth's atmosphere and ionosphere are assumed spherically symmetric, Bouguer's rule can be used to describe the refraction of a ray path in terms of a constant impact parameter (e.g., Budden, 1985),
where a is the impact parameter; r is the radial distance from the center of the curvature of the refracted ray to any point of the ray path; n is the refractive index (at radial distance r), which is related to refractivity N via $n=\mathrm{1}+{\mathrm{10}}^{\mathrm{6}}N$; and Φ is the local angle between the radial position vector and the ray direction at any point of the ray.
Equation (1) implies that, at each point along the ray path, the impact parameter a is equal to its initial value at the GPS transmitter position in the case of spherical symmetry, which leads to
where index i counts the (numerical 3D ray tracer) points along the ray path starting at the GPS transmitter position ${\overrightarrow{\mathit{R}}}_{\mathrm{G}}$ and ending at LEO; ${\overrightarrow{\mathit{R}}}_{i}$ and ${\widehat{\mathit{r}}}_{i}$ are the radial position vector and unit vector along the ray direction at point i, respectively; and Φ_{G} is the local angle between position vector and (initial) ray direction at the GPS transmitter, where we can furthermore assume that the refractivity is zero.
In reality nonsphericalsymmetry conditions of appreciable size will frequently occur, in particular between the ionospheric signal propagation inbound from the GPS and (after propagating through the atmosphere at tangent heights below 80 km) the one outbound to LEO (cf. Fig. 1); see, for example, the RO events discussed by Liu et al. (2013). In order to inspect the impact parameter changes along the ray path in these cases where a_{i} computed according to the second righthandside term of Eq. (2) will vary along the ray path, we cocompute the delta impact parameter Δa_{i} as the difference of the impact parameter at points i of the ray path and the impact parameter a_{G} at the GPS location:
Complementary to the geometrical parameters available from the ray tracing, the refractivity N comprises atmospheric and ionospheric terms, which are formulated based on standard equations as (e.g., Liu et al., 2015; Eqs. 1 and 4 therein)
where C_{atm}= 77.60 K hPa^{−1} and C_{ion}= 40.31 × 10^{6} m^{3} s^{−2} are the classical refractivity coefficients, p [hPa] and T [K] are atmospheric pressure and temperature (modeled by MSIS90), N_{e} [m^{−3}] is the ionospheric electron density (modeled by NeUoG), and f [Hz] is the GPS signal frequency (f_{L1}= 1.57542 GHz; f_{L2}= 1.22760 GHz).
In addition, the accumulated bending angle α_{i}, which accrues from the GPS position to any point i along the ray path, can be readily computed as the angle between the initial ray direction (unit vector ${\widehat{\mathit{r}}}_{\mathrm{G}})$ and the ray direction at point i (unit vector ${\widehat{\mathit{r}}}_{i})$:
The total bending angle along the entire ray is hence obtained by finally computing the angle between initial direction at GPS and terminal direction at LEO (subscript L):
Furthermore and importantly, the accumulated bending angle RIE, δα_{RIE}, can be estimated (after linearly interpolating in the alongray distance coordinate to the ray points i of the reference bending angle obtained without the ionosphere) by subtracting the reference bending angle α_{r} from the ionospherecorrected bending angle α_{c} (with the latter obtained by the standard dualfrequency correction of bending angles; e.g., Liu et al., 2015; Eq. 3 therein):
As a complement to these alongray accumulated quantities, the local bending angles and bending angle RIEs caused by individual ray tracer steps can also be readily cocomputed, by differencing the values between adjacent points i+1 and i:
Figures 4 to 7 sequentially illustrate for the three representative RO events (Occ.25, Occ.20, Occ.530) the alongray behavior of the key variables atmospheric and ionospheric refractivity (Eq. 4), L1 and L2 delta impact parameter (Eq. 3), L1 and L2 accumulated bending angle (Eq. 5), and ionospherecorrected Lc bending angle and bending angle RIE (Eq. 7). This is done in the form of imaging these variables for the three RO events in the alongray distance vs. impact height coordinate system (panels a and b of Figs. 4–7; ±3500 km alongray distance about ray tangent points; impact height range: 20 to 80 km) and in the form of depicting the alongray behavior of the two exceptional events along representative impact heights (80, 50, 30 km; panels c–f in Fig. 4 and panels c–d in Figs. 5–7).
In order to enable close inspection of the critical role of ionospheric symmetries, each of the Figs. 4–7 directly intercompares the nonspherical and spherical symmetry conditions. In terms of solar activity only the results for the medium solar activity level (F10.7 = 140) are illustrated, since we found that the influence of solar activity (which mainly drives the ionization level in the NeUoG model) is primarily to steer the magnitude of the effects (see Liu et al., 2013, 2015). The typical alongray characteristics are therefore reasonably well represented by just illustrating the mediumsolaractivity case.
Figure 4 shows the atmospheric and ionospheric refractivities and underpins that the alongray differences of inbound ionosphere (from the GPS) and outbound ionosphere (towards LEO) refractivities can be substantial. For example, in the case of the Occ.20 event these refractivities differ by more than a factor of 2 near the ionospheric F layer maximum, where the refractivities are largest (e.g., L2 inbound refractivity near 10 NU; L2 outbound refractivity reaching more than 20 NU). As expected, the atmospheric refractivity starts to exceed 1 NU only below about 35 km, and of course it exhibits no frequency dependence. It is thus essential to have a reliable firstorder and higherorder ionospheric correction to strongly mitigate the ionospheric effects that appear prominent down to the lower stratosphere.
Figure 5 shows the L1 and L2 delta impact parameters, which first of all verifies the reliability of the numerical ray tracing estimates, since the spherically symmetric ionosphere conditions indeed lead to alongray impact parameter changes of within 1 m only. This confirms that under such spherical symmetry conditions the bending angle retrievals (e.g., Schwarz et al., 2017) will be highly accurate, including for the impact height that is decisive for enabling accurate vertical geolocation (ScherllinPirscher et al., 2017). Under nonsphericalsymmetry ionosphere conditions, the Occ.20 event with the largest asymmetry of the example events shows that alongray L1 and L2 impact parameter variations of more than 10 to 20 m are possible and are generally found to be negative (relative to the initial impact parameter at the GPS transmitter). Lc bending angle retrievals with their intrinsic spherical symmetry assumption should thus receive higherorder ionospheric correction to mitigate such possible impacts.
Figure 6 depicts the accumulated L1 and L2 bending angles, which highlight the significant alongray modulations that the bending angle receives due to the ionospheric influences relative to the atmospheric bending angle, in particular above about 35 km in the upper stratosphere and mesosphere, where the neutral atmospheric bending angle is rather small. Below 35 km the dominating influence of the atmosphere in the vicinity of the tangent point location becomes prominently visible, in line with the exponential increase of atmospheric refractivity (Fig. 4) and hence atmospheric bending down into the lower stratosphere. Nevertheless even at 30 km the ionospheric contribution is still visible, which underscores that an accurate ionospheric correction with minimized residual error will be vital.
Figure 7a shows the ionospherecorrected Lc bending angle and indicates, compared to Figs. 6a and b, that the standard linear dualfrequency correction of bending angles basically does a very effective job in eliminating the ionospheric bending angle contributions. The Lc bending angle images look visually very clean and are highly dominated by just the atmospheric accumulated bending angle accruing at all heights around the tangent point location. Directly inspecting the bending angle RIE, finally, shows that the alongray behavior and accumulated magnitude of the higherorder RIE left by the linear correction significantly depend in particular on asymmetry conditions. Technical ray tracer effects are also visible as intermittent spiky behavior, since the RIEs are at the submicroradian to microradian magnitude level only, which is a challenge for the ionospheric model smoothness as discussed in Sect. 2.1.
The Occ.530 event under nonsphericalsymmetry appears to accumulate the highest RIEs of near 2 to 4 µrad at LEO, while the spherical symmetry cases both accumulate RIEs up to around 0.5 to 1 µrad or less only. This is in line with Fig. 2, which shows for the majority of the 26 exceptional events the dominance of asymmetry and 3D ray tracing effects in driving the RIE magnitude. Also, as shown by Fig. 7c and d (and found for other RO events inspected but not separately shown), the mesospheric RIEs above about 50 km generally appear to be higher than the upperstratospheric ones from 50 km downwards. This is in line with findings of Syndergaard (2000) and likely driven by the increased closeness of the tangent point height to the ionospheric E layer peaking near 105 km, which makes the Lc bending angle more vulnerable to higherorder RIEs.
Figure 7c and d (and alongray results for further exceptional events not separately shown) also clearly indicate the mixingin of technical ray tracer effects in our simulations. These render it more difficult to rigorously quantify the (smooth) physical RIE effects from the largescale ionospheric model structure since, despite the reasonable smoothing applied, the spiky components may somewhat perturb the smooth accumulated results as well. In the future we will therefore aim to further improve the simulation setup to fully isolate the technical from the physicsbased propagation effects. For now we have found clear evidence, nevertheless, that currently both technical effects and cases with physically high RIEs from ionospheric asymmetries play important roles in explaining the anomalous behavior of the exceptional RO events.
Previous theoretical and simulation studies as well as empirical studies that we surveyed in the Introduction have characterized and quantified higherorder residual ionospheric errors (RIEs) in bending angles by analyses of individual events as well as ensembles of events. The statistical results showed that the mean bending angle RIE biases are predominantly negative, typically at the 0.03 to 0.1 µrad level, and these biases may lead to systematic errors in stratospheric climatologies built from retrieved profiles. The RIE standard deviations are typically at the 0.1 to 0.5 µrad level, and they have a clear tendency to increase with increasing solar activity, i.e., with increasing ionization level (electron density) in the ionosphere.
In our previous Liu et al. (2015) study we had contributed to these findings but were left with 26 exceptional RO events with very high RIEs, at the 1 to 10 µrad standard deviation level, in the context of about 700 standard events with lownoise RIEs within 0.5 µrad standard deviation. In this study we therefore placed focus on these 26 exceptional events and, by way of detailed alongray analyses of ray tracing simulations over the stratosphere and mesosphere, inspected the causes of anomalously high RIEs. The goal at the same time was to deepen quantitative insight into how RIEs accumulate during signal propagation, along with accumulation of the total atmospheric bending angles that are the desired RO observables.
From the results of these analyses we conclude with the following main findings on the causes of the exceptional RO events:

Strengthening previous results by Mannucci et al. (2010, 2011), we find that asymmetric ionospheric conditions play an important role for anomalously high RIEs, more than ionization levels driven by solar activity and possible geographic location dependencies that seemed to be present from salient geographic clustering of the majority of exceptional RO events in two regions (European–Asian region and Indian Ocean region).

The fact that no obvious physicsbased explanation was found for the geographic clustering and the intermittent spiky behavior found in simulated RIEs indicates that a portion of the anomalous RIEs of the exceptional RO events were caused by the technical challenge of ray tracing at millimetric excessphase accuracy through the 3D ionospheric model NeUoG, which is not perfectly smooth everywhere in its electron density field derivatives.

The detailed alongray analyses of atmospheric and ionospheric refractivities, impact parameter changes, bending angles, and RIEs also revealed that alongray L1 and L2 impact parameter variations of more than 10 to 20 m are possible due to ionospheric asymmetries and are generally found to be negative (relative to the initial impact parameter at the GPS transmitter). Standard bending angle retrievals with their intrinsic spherical symmetry assumption should thus receive higherorder ionospheric correction to mitigate such impacts.

The mesospheric RIEs above about 50 km generally appear to be higher than the upperstratospheric ones from 50 km downwards. This is in line with findings of Syndergaard (2000) and likely driven by the increased closeness of the tangent point height to the ionospheric E layer peaking near 105 km, which makes the standard ionospherecorrected bending angles more vulnerable to higherorder RIEs.
Overall this study of exceptional RO events with anomalous RIEs in our endtoend simulations indicated that the main causes are a combination of physicsbased effects, in particular ionospheric asymmetries, and of technical ray tracer effects due to occasionally imperfect smoothness of modeling ionospheric refractivity field derivatives. This makes it more difficult to rigorously quantify the physicsbased RIE effects from the largescale ionospheric model structure since the intermittent spiky nature of the technical effects may somewhat perturb the smooth accumulated results as well.
In the future we will therefore aim to further improve our alongray simulation and analysis system to fully isolate the technical from the physicsbased propagation effects. For now we have found clear evidence, nevertheless, that currently both technical effects and cases with physically high RIEs from ionospheric asymmetries play major roles in explaining the anomalous behavior of the exceptional RO events. The detailed alongray modeling system will also be valuable beyond this work for additional GNSS RO signal propagation studies.
The ECMWF (Reading, UK) is thanked for access to their archived analysis and forecast data (more information available at http://www.ecmwf.int/en/ forecasts/datasets). The software code used for this study does not belong to the public domain and cannot be distributed. To access the relevant result files of this study, please contact the corresponding author.
The authors declare that they have no conflict of interest.
This article is part of the special issue “Observing Atmosphere and Climate with Occultation Techniques – Results from the OPACIROWG 2016 Workshop”. It is a result of the International Workshop on Occultations for Probing Atmosphere and Climate, Leibnitz, Austria, 8–14 September 2016.
This research was partially supported by the National Natural Science
Foundation of China (grant nos. 41405039, 41775034, 41405040, 41505030,
41606206, and 41730109) and the FengYun3 (FY3) Global Navigation Satellite
System Occultation Sounder (GNOS) development and manufacture project led by
NSSC, CAS. The research at WEGC/University of Graz was supported by the European
Space Agency (ESA) projects OPSGRAS and MMValRO and the Austrian Research
Promotion Agency (FFG) project OPSCLIMPROP (ASAP9 project no. 840070). We
acknowledge Johannes Fritzer (WEGC) for his support in EGOPS software developments
valuable to this study. The research at SPACE/RMIT University was supported
by the Australian Research Council (ARC) (LP0883288), the Australian
Antarctic Division (project no. 4159), and the CAS/SAFEA International
Partnership Program for Creative Research Teams (grant no. KZZDEWTZ05).
The support from the Jiangsu dual creative talents and Jiangsu dual creative
team program projects awarded to CUMT in 2017 is acknowledged.
Edited by: Sean Healy
Reviewed by: two anonymous referees
Angerer, B., Ladstädter, F., ScherllinPirscher, B., Schwärz, M., Steiner, A. K., Foelsche, U., and Kirchengast, G.: Quality aspects of the Wegener Center multisatellite GPS radio occultation record OPSv5.6, Atmos. Meas. Tech., 10, 4845–4863, https://doi.org/10.5194/amt1048452017, 2017.
Angling, M. J., Elvidge, S., and Healy, S. B.: Improved model for correcting the ionospheric impact on bending angle in radio occultation measurements, Atmos. Meas. Tech., 11, 2213–2224, https://doi.org/10.5194/amt1122132018, 2018.
Anthes, R. A.: Exploring Earth's atmosphere with radio occultation: contributions to weather, climate and space weather, Atmos. Meas. Tech., 4, 1077–1103, https://doi.org/10.5194/amt410772011, 2011.
Bassiri, S. and Hajj, G. A.: Higherorder ionospheric effects on the GPS observables and means of modeling them, Manuscr. Geodaet., 18, 280–289, 1993.
Budden, K. G.: The propagation of radio waves: the theory of radio waves of low power in the ionosphere and magnetosphere, Cambridge University Press, Cambridge, 1985.
Coleman, C. J. and Forte, B.: On the residual ionospheric error in radio occultation measurements, Radio Sci., 52, 918–937, https://doi.org/10.1002/2016RS006239, 2017.
Danzer, J., ScherllinPirscher, B., and Foelsche, U.: Systematic residual ionospheric errors in radio occultation data and a potential way to minimize them, Atmos. Meas. Tech., 6, 2169–2179, https://doi.org/10.5194/amt621692013, 2013.
Danzer, J., Healy, S. B., and Culverwell, I. D.: A simulation study with a new residual ionospheric error model for GPS radio occultation climatologies, Atmos. Meas. Tech., 8, 3395–3404, https://doi.org/10.5194/amt833952015, 2015.
Edwards, P. G. and Pawlak, D.: Metop: The space segment for Eumetsat's Polar System, ESA Bull., 102, 6–18, 2000.
Foelsche, U., ScherllinPirscher, B., Ladstädter, F., Steiner, A. K., and Kirchengast, G.: Refractivity and temperature climate records from multiple radio occultation satellites consistent within 0.05 %, Atmos. Meas. Tech., 4, 2007–2018, https://doi.org/10.5194/amt420072011, 2011.
Fritzer, J., Kirchengast, G., and Pock, M.: EndtoEnd Generic Occultation Performance Simulation and Processing System version 5.6 (EGOPS 5.6) Software User Manual, Tech. Rep. ESAESTEC WEGCEGOPS2013TR01, Wegener Center and Inst. for Geophys., Astrophys., and Meteorol., Univ. of Graz, Austria, 2013.
Gobiet, A. and Kirchengast, G.: Advancements of Global Navigation Satellite System radio occultation retrieval in the upper stratosphere for optimal climate monitoring utility, J. Geophys. Res., 109, D24110, https://doi.org/10.1029/2004jd005117, 2004.
Gobiet, A., Kirchengast, G., Manney, G. L., Borsche, M., Retscher, C., and Stiller, G.: Retrieval of temperature profiles from CHAMP for climate monitoring: intercomparison with Envisat MIPAS and GOMOS and different atmospheric analyses, Atmos. Chem. Phys., 7, 3519–3536, https://doi.org/10.5194/acp735192007, 2007.
Gorbunov, M. E.: Ionospheric correction and statistical optimization of radio occultation data, Radio Sci., 37, 1084, https://doi.org/10.1029/2000rs002370, 2002.
Hajj, G. A., Kursinski, E. R., Romans, L. J., Bertiger, W. I., and Leroy, S. S.: A technical description of atmospheric sounding by GPS occultation, J. Atmos. Sol.Terr. Phy., 64, 451–469, https://doi.org/10.1016/S13646826(01)001146, 2002.
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.
Hedin, A. E.: Extension of the MSIS thermosphere model into the middle and lower atmosphere, J. Geophys. Res., 96, 1159–1172, https://doi.org/10.1029/90JA02125, 1991.
Ho, S.P., Hunt, D., Steiner, A. K., Mannucci, A. J., Kirchengast, G., Gleisner, H., Heise, S., von Engeln, A., Marquardt, C., Sokolovskiy, S., Schreiner, W., ScherllinPirscher, B., Ao, C., Wickert, J., Syndergaard, S., Lauritsen, K. B., Leroy, S., Kursinski, E. R., Kuo, YH., Foelsche, U., Schmidt, T., and Gorbunov, M.: Reproducibility of GPS radio occultation data for climate monitoring: Profiletoprofile intercomparison of CHAMP climate records 2002 to 2008 from six data centers, J. Geophys. Res., 117, D18111, https://doi.org/10.1029/2012JD017665, 2012.
Hoque, M. M. and Jakowski, N.: Higher order ionospheric propagation effects on GPS radio occultation signals, Adv. Space Res., 46, 162–173, https://doi.org/10.1016/j.asr.2010.02.013, 2010.
Hoque, M. M. and Jakowski, N.: Ionospheric bending correction for GNSS radio occultation signals, Radio Sci., 46, RS0D06, https://doi.org/10.1029/2010rs004583, 2011.
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., 102, 23429–23465, https://doi.org/10.1029/97jd01569, 1997.
Lackner, B. C., Steiner, A. K., Hegerl, G. C., and Kirchengast, G.: Atmospheric climate change detection by radio occultation data using a fingerprinting method, J. Climate, 24, 5275–5291, https://doi.org/10.1175/2011JCLI3966.1, 2011.
Ladreiter, H. P. and Kirchengast, G.: GPS/GLONASS sensing of the neutral atmosphere: Modelindependent correction of ionospheric influences, Radio Sci., 31, 877–891, https://doi.org/10.1029/96rs01094, 1996.
Leitinger, R. and Kirchengast, G.: Easy to use global and regional ionospheric models – A report on approaches used in Graz, Acta Geod. Geophys. Hu., 32, 329–342, https://doi.org/10.1007/BF03325504, 1997.
Liu, C. L., Kirchengast, G., Zhang, K. F., Norman, R., Li, Y., Zhang, S. C., Carter, B., Fritzer, J., Schwaerz, M., Choy, S. L., Wu, S. Q., and Tan, Z. X.: Characterisation of residual ionospheric errors in bending angles using GNSS RO endtoend simulations, Adv. Space Res., 52, 821–836, https://doi.org/10.1016/j.asr.2013.05.021, 2013.
Liu, C. L., Kirchengast, G., Zhang, K., Norman, R., Li, Y., Zhang, S. C., Fritzer, J., Schwaerz, M., Wu, S. Q., and Tan, Z. X.: Quantifying residual ionospheric errors in GNSS radio occultation bending angles based on ensembles of profiles from endtoend simulations, Atmos. Meas. Tech., 8, 2999–3019, https://doi.org/10.5194/amt829992015, 2015.
Loiselet, M., Stricker, N., Menard, Y., and Luntama, J.P.: GRAS – Metop's GPSbased atmospheric sounder, ESA Bull., 102, 38–44, 2000.
Luntama, J.P., Kirchengast, G., Borsche, M., Foelsche, U., Steiner, A., Healy, S., von Engeln, A., O'Clerigh, E., and Marquardt, C.: Prospects of the EPS GRAS mission for operational atmospheric applications, B. Am. Meteorol. Soc., 89, 1863–1875, https://doi.org/10.1175/2008BAMS2399.1, 2008.
Mannucci, A. J., Ao, C. O., Pi, X., and Iijima, B. A.: Impact of the ionosphere on GNSS radio occultation retrievals, Presentation at the OPACGRASSAFIROWG International Workshop September 6–11, 2010, Graz Austria, available at: http://wegcwww.unigraz.at/opac2010/pdf_presentation/opac_2010_mannucci_anthony_presentation71.pdf (last access: 23 March 2018), 2010.
Mannucci, A. J., Ao, C. O., Pi, X., and Iijima, B. A.: The impact of large scale ionospheric structure on radio occultation retrievals, Atmos. Meas. Tech., 4, 2837–2850, https://doi.org/10.5194/amt428372011, 2011.
Melbourne, W. G., Davis, E. S., Duncan, C. B., Hajj, G. A., Hardy, K. R., Kursinski, E. R., Meehan, T. K., Yong, L. E., and Yunck, T. P.: The application of spaceborne GPS to atmospheric limb sounding and global change monitoring, JPL Publication 9418, Jet Propulsion Lab., Cal. Tech., Pasadena, CA, 1994.
ScherllinPirscher, B., Kirchengast, G., Steiner, A. K., Kuo, Y.H., and Foelsche, U.: Quantifying uncertainty in climatological fields from GPS radio occultation: an empiricalanalytical error model, Atmos. Meas. Tech., 4, 2019–2034, https://doi.org/10.5194/amt420192011, 2011.
ScherllinPirscher, B., Steiner, A. K., Kirchengast, G., Schwaerz, M., and Leroy, S. S.: The power of vertical geolocation of atmospheric profiles from GNSS radio occultation, J. Geophys. Res. Atmos., 122, 1595–1616, https://doi.org/10.1002/2016JD025902, 2017.
Schwarz, J., Kirchengast, G., and Schwaerz, M.: Integrating uncertainty propagation in GNSS radio occultation retrieval: from excess phase to atmospheric bending angle profiles, Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt2017159, in review, 2017.
Steiner, A. K. and Kirchengast, G.: Error analysis for GNSS radio occultation data based on ensembles of profiles from endtoend simulations, J. Geophys. Res., 110, 1–21, https://doi.org/10.1029/2004JD005251, 2005.
Steiner, A. K., Kirchengast, G., Lackner, B. C., Pirscher, B., Borsche, M., and Foelsche, U.: Atmospheric temperature change detection with GPS radio occultation 1995 to 2008, Geophys. Res. Lett., 36, L18702, https://doi.org/10.1029/2009GL039777, 2009.
Steiner, A. K., Lackner, B. C., Ladstädter, F., ScherllinPirscher, B., Foelsche, U., and Kirchengast, G.: GPS radio occultation for climate monitoring and change detection, Radio Sci., 46, RS0D24, https://doi.org/10.1029/2010RS004614, 2011.
Steiner, A. K., Hunt, D., Ho, S.P., Kirchengast, G., Mannucci, A. J., ScherllinPirscher, B., Gleisner, H., von Engeln, A., Schmidt, T., Ao, C., Leroy, S. S., Kursinski, E. R., Foelsche, U., Gorbunov, M., Heise, S., Kuo, Y.H., Lauritsen, K. B., Marquardt, C., Rocken, C., Schreiner, W., Sokolovskiy, S., Syndergaard, S., and Wickert, J.: Quantification of structural uncertainty in climate data records from GPS radio occultation, Atmos. Chem. Phys., 13, 1469–1484, https://doi.org/10.5194/acp1314692013, 2013.
Syndergaard, S.: Modeling the impact of the Earth's oblateness on the retrieval of temperature and pressure profiles from limb sounding, J. Atmos. Sol.Terr. Phys., 60, 171–180, https://doi.org/10.1016/S13646826(97)000564, 1998.
Syndergaard, S.: Retrieval analysis and methodologies in atmospheric limb sounding using the GNSS radio occultation technique, DMI Sci. Rep. 996, Danish Meteorol. Inst., Copenhagen, Denmark, 1999.
Syndergaard, S.: On the ionosphere calibration in GPS radio occultation measurements, Radio Sci., 35, 865–883, https://doi.org/10.1029/1999rs002199, 2000.
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, Phys. Atmos. Ocean, 29, 602–609, 1994.