Quantifying residual ionospheric errors in GNSS radio occultation bending angles based on ensembles of profiles from end-to-end simulations

The radio occultation (RO) technique using signals from the Global Navigation Satellite System (GNSS), in particular from the Global Positioning System (GPS) so far, is currently widely used to observe the atmosphere for applications such as numerical weather prediction and global climate monitoring. The ionosphere is a major error source in RO measurements at stratospheric altitudes, and a linear ionospheric correction of dual-frequency RO bending angles is commonly used to remove the first-order ionospheric effect. However, the residual ionospheric error (RIE) can still be significant so that it needs to be further mitigated for high-accuracy applications, especially above about 30 km altitude where the RIE is most relevant compared to the magnitude of the neutral atmospheric bending angle. Quantification and careful analyses for better understanding of the RIE is therefore important for enabling benchmark-quality stratospheric RO retrievals. Here we present such an analysis of bending angle RIEs covering the stratosphere and mesosphere, using quasi-realistic end-to-end simulations for a full-day ensemble of RO events. Based on the ensemble simulations we assessed the variation of bending angle RIEs, both biases and standard deviations, with solar activity, latitudinal region and with or without the assumption of ionospheric spherical symmetry and co-existing observing system errors. We find that the bending angle RIE biases in the upper stratosphere and mesosphere, and in all latitudinal zones from low to high latitudes, have a clear negative tendency and a magnitude increasing with solar activity, which is in line with recent empirical studies based on real RO data although we find smaller bias magnitudes, deserving further study in the future. The maximum RIE biases are found at low latitudes during daytime, where they amount to within −0.03 to −0.05 μrad, the smallest at high latitudes (0 to −0.01 μrad; quiet space weather and winter conditions). Ionospheric spherical symmetry or asymmetries about the RO event location have only a minor influence on RIE biases. The RIE standard deviations are markedly increased both by ionospheric asymmetries and increasing solar activity and amount to about 0.3 to 0.7 μrad in the upper stratosphere and mesosphere. Taking also into account the realistic observation errors of a modern RO receiving system, amounting globally to about 0.4 μrad (unbiased; standard deviation), shows that the random RIEs are typically comparable to the total observing system error. The results help to inform future RIE mitigation schemes that will improve upon the use of the linear ionospheric correction of bending angles and also provide explicit uncertainty estimates. Published by Copernicus Publications on behalf of the European Geosciences Union. 3000 C. L. Liu et al.: Quantifying residual ionospheric errors in bending angles

the upper stratosphere and mesosphere, and in all latitudinal zones from low to high latitudes, have a clear negative tendency and a magnitude increasing with solar activity, which is in line with recent empirical studies based on real RO data although we find smaller bias magnitudes, deserving further study in the future.The maximum RIE biases are found at low latitudes during daytime, where they amount to within −0.03 to −0.05 µrad, the smallest at high latitudes (0 to −0.01 µrad; quiet space weather and winter conditions).Ionospheric spherical symmetry or asymmetries about the RO event location have only a minor influence on RIE biases.The RIE standard deviations are markedly increased both by ionospheric asymmetries and increasing solar activity and amount to about 0.3 to 0.7 µrad in the upper stratosphere and mesosphere.Taking also into account the realistic observation errors of a modern RO receiving system, amounting globally to about 0.4 µrad (unbiased; standard deviation), shows that the random RIEs are typically comparable to the total observing system error.The results help to inform future RIE mitigation schemes that will improve upon the use of the linear ionospheric correction of bending angles and also provide explicit uncertainty estimates.

Introduction
Detection of global climate change is a significant challenge in atmospheric sciences (Steiner et al., 2011) due to the extreme complexity and dynamics of the Earth's atmospheric system (Zhang et al., 2011) and to the stringent climate monitoring principles such as reproducibility, homogeneity, longterm stability, high accuracy, high spatial and temporal resolution and global coverage.In addition to these principles, the Global Climate Observing System 2010 (GCOS, 2010) defined observation requirements for essential climate variables, such as upper-air tropospheric and stratospheric temperature.The requirements for the precision and resolution of temperature profiles are less than 0.5 K root-mean-square value, 500 and 0.5 km horizontal and vertical resolutions, respectively, in the upper troposphere and 1.5 km vertical resolution in the lower stratosphere (Immler et al., 2010;Steiner et al., 2011).
Current atmospheric observation techniques such as radiosondes and weather satellite radiometers can hardly meet these requirements.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 remote sensing technique.It can deliver data traceable to the international standard of time (the SI second) and has the potential for monitoring decadal-scale climate change (Steiner et al., 2009;Scherllin-Pirscher et al., 2011b;Lackner et al., 2011) due to its unique characteristics such as high vertical resolution, high accuracy and long-term stability of its observations, as well as self-calibration capability and global coverage (Gobiet and Kirchengast, 2004;Steiner et al., 2011).Figure 1 illustrates the GPS-to-LEO occultation geometry.
Detailed analyses of GNSS RO errors have been conducted by many scientists (Kursinski et al., 1997;Rieder and Kirchengast, 2001;Steiner and Kirchengast, 2005;Scherllin-Pirscher et al., 2011a, b).These errors mainly include the satellites' orbital error, clock biases, systematic hardware delay, antenna phase center variation, cycle slips, ionospheric refraction, atmospheric multipath and scintillations.These studies demonstrated that GNSS RO observations have the best quality in the upper troposphere and lower stratosphere (UTLS, defined as the 5-35 km altitude range).
The accuracy of atmospheric variables retrieved from GPS/MET observations was found to be 0.4 % in refractivity and 1 K in temperature, at 0.5-1 km vertical resolution in the UTLS (Kursinski et al., 1996;Rocken et al., 1997;Steiner et al., 1999).The climatological analyses of CHAMP data revealed an accuracy of 0.4 K in the global-mean temperature in the UTLS and a 1 K standard deviation at an altitude of 10 km that increased to 2 K at 30 km (Wickert et al., 2004).A comparison between COSMIC RO results and co-located radiosonde data showed high agreement in temperature profiles, with less than 0.5 K mean differences and less than 2.0 K standard deviations of all the bins (He et al., Figure 1.Radio occultation geometry illustrating the separated L1 and L2 signal paths and the ionosphere-corrected ray path L c ; α c is the ionosphere-corrected bending angle, a is the impact parameter and r is the radius from the Earth's center of curvature to the tangent point of the GPS-LEO signal path. 2009).Foelsche et al. (2009) showed that the monthly means of CHAMP, GRACE-A and COSMIC global-average climatology agreed well, with a < 0.05 % discrepancy in refractivity and < 0.05 % in dry temperature for almost all months in the UTLS (Foelsche et al., 2011).Ho et al. (2012) and Steiner et al. (2013) demonstrate the very low structural uncertainty of the data.
However, these studies also showed that the errors in the RO retrievals above about 35 km were significant (Rieder and Kirchengast, 2001;Gobiet et al., 2007;Scherllin-Pirscher et al., 2011b).This is because from the stratosphere to the mesosphere, ionospheric error influences become more and more dominant (e.g., Mannucci et al., 2011;Liu et al., 2013).A correction for the ionospheric effects in RO is usually implemented using a dual-frequency linear combination of bending angles (Vorob'ev and Krasil'nikova, 1994;Syndergaard, 2000), and the remaining error is the so-called residual ionospheric error (RIE).The RIE is not negligible in the upper stratosphere and mesosphere (USMS) in climate applications (Syndergaard, 2000;Mannucci et al., 2011;Danzer et al., 2013;Liu et al., 2013).What is worse is that the RIE can propagate downwards into the UTLS atmospheric retrievals through the Abel integral and the hydrostatic integral (e.g., Kursinski et al., 1997;Steiner and Kirchengast, 2005).
Due to the fact that the RIE in atmospheric profiles at high altitudes is large, a climatological model such as MSIS-90 is often used to obtain atmospheric variable values for high altitudes instead of using RO atmospheric retrievals, and the model-derived values are so-called background information.This background information is commonly used for the initialization of high-altitude atmospheric profiles.For the part in mid-to-high altitudes, the background information, more specifically, the background bending angles, together with RO-derived bending angles, often called observed bending angles, is used in a statistical optimization approach to obtain optimal bending angle profiles (Sokolovskiy and Hunt, 1996;Kursinski et al., 1997;Hocke, 1997;Steiner et al., 1999;Healy, 2001;Gorbunov, 2002;Gobiet and Kirchengast, 2004;Li et al., 2013).
Statistical optimization does not improve the quality of the observed bending angle profiles (Gobiet and Kirchengast, 2004) but instead just helps to better initialize the Abel integration.Thus it is necessary to develop a better ionospheric calibration scheme for further improvement of the RO atmospheric bending angle retrievals.That is, in order to improve the quality of RO retrievals in the USMS for extending RO observations' climatological utility up to or exceeding the stratopause, the characterization of the RIE is significant for effective mitigation of its effect on RO retrievals.
In this study, based on the end-to-end simulation approach, an ensemble simulation using a realistic transmitter-receiver geometry is conducted to characterize and quantify RIEs in daily-global-mean bending angle profiles.The 3-D NeUoG ionospheric model (Leitinger and Kirchengast, 1997) and the MSIS-90 neutral atmospheric model were used in the simulations.In Sect.2, our simulation strategy for analyzing the bending angle RIEs will be elaborated, followed by an ensemble simulation design and results of analyses in Sect.3. A summary and conclusions are finally given in Sect. 4.

Bending angle RIEs and simulation method
The Earth's ionosphere is a mixed neutral-and-ionized gas consisting mainly of free electrons, ions, neutral atoms and molecules, located at the altitude range of around 80-1000 km.Since the ionosphere is a dispersive medium, the magnitudes of GNSS signal carrier phase delays and bending angles are related to their frequency.The first-order ionospheric effect can be largely mitigated by a dual-frequency linear combination of GPS signal observations.The residual ionospheric errors in GNSS RO retrievals contain not only the residual first-order effect but also high-order effects, caused by the uneven distribution and anisotropy of the ionospheric plasma, respectively.
Early studies of the RIE for GNSS geodetic applications have been conducted by several researchers (Brunner and Gu, 1991;Bassiri and Hajj, 1993;Hoque and Jakowski, 2007).They found that the RIEs in the context of space-toground GNSS positioning uses are mainly contributed by the high-order effect and the bending effect of the signal carriers.Since the emergence of the GNSS RO concept, the RIE effects on the GNSS RO retrievals have also been investigated by several scientists.For example, Ladreiter and Kirchengast (1996) took into account the splitting effect of GPS dual-frequency signals and suggested a model-independent ionospheric calibration approach similar to Vorob'ev and Krasil'nikova (1994).Syndergaard (2000) performed a detailed theoretical analysis of the dual-frequency signals' bending and splitting effects on the RIEs and proposed an improved phase path correction method.His study showed that the first-order dispersion residual is dominant in the RIEs.Hoque and Jakowski (2010) used the ray tracing method to study the effects of the high-order ionospheric propagation on GPS radio occultation signals.They found that the estimated maximum separation of the dual-frequency ray paths reached 1 km, and the maximum excess phase was about 2.7 m in which the second-and third-order ionospheric effects were about 13 and 2.7 cm, respectively.Mannucci et al. (2011) conducted a simulation study to assess the magnitude of the RIEs under ionospheric storm conditions by a study of the propagation of GPS signals in an occulting geometry.They concluded that RO retrievals suffer from ionospheric storms dramatically, especially at the heights above 25 km.
Recently Liu et al. (2013) performed an initial study on quantifying bending angle RIEs based on end-to-end simulations, complemented by Liu et al. (2014) who looked at the effects in subsequently retrieved temperature profiles.This work found, based on a limited set of individual occultation events, that the RIEs can significantly affect bending angles and retrieved temperatures, in particular when ionospheric disturbances occur such as during periods of active space weather (bending angle RIEs can exceed 0.5 µrad and temperature errors 1 K in the upper stratosphere).These initial results provided the encouragement and formed the basis for the present, much more advanced, study on bending angle RIEs based on much larger ensembles.

Ionospheric correction and bending angle RIE
The refractive index of the ionosphere and ionospheric radio wave propagation in detail can be found in the literatures (e.g., Budden, 1985;Davies, 1990;Brunner and Gu, 1991;Bassiri and Hajj, 1993;Ladreiter and Kirchengast, 1996).For GPS signals the series expansion of the ionospheric refractive index can be approximated by where the two constants C = e 2 /(8π 2 mε 0 ) ∼ = 40.308m 3 s −2 and K = e 3 /(16π 3 m 2 ε 0 ) ∼ = 1.1283 × 10 12 m 3 /(Ts 3 ), N e is the electron density, B is the geomagnetic field strength, f is the radio wave frequency, e is the elementary charge, m is the electron mass, ε 0 is the permittivity of vacuum and θ is the angle between the magnetic field vector and the wave propagation direction.For the GPS signal frequencies L1 and L2 (f 1 = 1.57542GHz and f 2 = 1.22760GHz), the magnitudes of the first-and second-order terms on the right-hand side of Eq. (1) are 10 −4 and 10 −7 , respectively.The second-order geomagnetic term can be generally neglected in GPS RO applications (Syndergaard, 2000;Hoque and Jakowski, 2010;Liu et al., 2013), although it is a main source of the RIE in space-to-ground GPS positioning applications (Hoque and Jakowski, 2007).Due to this, as well as for computational efficiency, this term was ignored in our global ensemble RO data simulation for this study.We suggest that further refined RIE studies should revisit the role of the geomagnetic term again, however, since it remains unclear whether small residual biases (possibly within 0.01 to 0.1 µrad) may arise under high ionization levels in some geographic regions despite the general smallness of the term.Such a quantification will be one useful element for helping reconcile a current discrepancy in bending angle RIE bias estimations from simulations like here and RIE estimations from real RO data such as by Danzer et al. (2013); see further comments in Sect. 4 below.
Electron density N e is the key physical quantity of the ionosphere due to its dominant effect on radio wave propagation.As GPS signal carriers are in the L band, the effect of free electrons on the refractive index is larger than that of the neutral gas per unit mass (Mannucci et al., 2011).The maximum daytime ionospheric refractive index occurs at a height of around 300 km and is comparable in magnitude to the atmospheric refractive index at a height of around 20-30 km.For high-accuracy atmospheric variable retrievals, this magnitude of ionospheric effect clearly needs to be corrected.The most common approach to correcting the ionospheric effect is to use a dual-frequency linear combination of GPS excess phase observations in GPS data processing as expressed by (Spilker and Jr, 1978): where L i (t)(i = 1, 2) is the optical phase path observations between the transmitter and the receiver, t is the sample time/epoch and L c (t) is the dual-frequency linear combination of the L1 and L2 phase observations.The combination of phase path observations was used in the early stage of GPS RO, as widely used in GPS geodetic applications.In ground-based GPS geodetic measurements, high-order ionospheric effects are dominant in the RIEs as mentioned above.However, in the RIEs of GPS RO observations the first-order residual is dominant due to the bending and ray path separation of the two frequency signals L1 and L2 (Syndergaard, 2000;Hoque and Jakowski, 2007, 2010, 2011), especially during the daytime and when a solar maximum period prevails (Syndergaard, 2000;Hoque and Jakowski, 2010;Mannucci et al., 2011).
To mitigate the separation effect of GPS carriers on the RIEs, Vorob'ev and Krasil'nikova (1994) proposed the dualfrequency linear combination of bending angles retrieved from the two frequencies at a common impact parameter, as expressed by where α 1 (a) and α 2 (a) are the two bending angles derived from the two frequency signals at the impact parameter a, and α c (a) is the dual-frequency linear combination of α 1 (a) and α 2 (a).
The combination of two bending angles is more effective than the phase combination expressed by Eq. (2) in the RO context (Hocke, 1997;Kursinski et al., 1997;Rocken et al., 1997;Gorbunov and Gurvich, 1998;Feng and Herman, 1999;Steiner et al., 1999), since it not only accounts for the separation of the dual-frequency carriers but also considers the fact that most of the total bending angle is accumulated near the ray perigee (Vorob'ev and Krasil'nikova, 1994;Ladreiter and Kirchengast, 1996).Due to its smaller ionospheric residual, the bending angle correction has become the most popular ionospheric correction approach nowadays in GPS RO data processing.Our simulation and quantification study therefore also focuses on the RIEs relative to this combination.
The bending angle RIE contaminates the accuracy of atmospheric profile retrievals when they are not corrected for.For high-accuracy meteorological monitoring and benchmark climate applications, more effective algorithms or approaches for mitigating the effect of the bending angle RIE are critical.Effective algorithms must be based on the characteristics of the actual errors.To investigate the characteristics of the global ensemble bending angle RIE, in this study, realistic simulations for bending angle RIEs using the ray tracing technique were conducted in which the 3-D NeUoG ionospheric model and the MSIS-90 neutral atmospheric model were used as the ionospheric and atmospheric background models.

RO end-to-end simulation tool
The End-to-end Generic Occultation Performance Simulation and Processing System (EGOPS) (Fritzer et al., 2011) was used as the simulation tool in this study, in which the whole process of simulating an RO event consists of the following five stages: (1) simulating satellite geometry, (2) modeling the neutral atmosphere and ionosphere, (3) simulating GPS dual-frequency signals' propagation through the atmosphere, (4) simulating the observation system and observations and (5) retrieving bending angles and other atmospheric profiles such as refractivity and temperature.

Atmospheric and ionospheric modeling
The focus of this study is to investigate the characteristics of bending angle RIEs in the USMS rather than the errors caused by the neutral atmosphere.Therefore, the MSIS-90 (Hedin, 1991) climatological atmospheric model was used in the ray tracing process due to its simplicity.Furthermore, to derive reliable bending angle RIEs, which are based on the dual-frequency linear combination of bending angles without the effects of the neutral atmosphere horizontal gradient and water vapor density ambiguity, the assumptions of local spherical symmetry and dry atmosphere were made in the simulation which are well justified for the focus heights of interest above 20 km.In this case, the neutral atmospheric refractivity, which depends only on atmospheric pressure P and temperature T , can be expressed as N = 77.6P/T . (4) The 3-D NeUoG ionospheric model (Leitinger and Kirchengast, 1997), as a function of location (latitude, longitude and height), month, universal time (UT) and solar activity (F10.7), was used for simulating the distribution of 3-D electron density.To investigate the effect of solar activity level on the characteristics and quantities of the bending angle RIEs, we simulated 3-D electron densities under low, medium and high solar activity levels by setting 70, 140 and 210 as the value of the solar radio flux index F10.7.Then the ionospheric refractivity was calculated by the approximated ionospheric refractivity formula (Eq.( 1) without the magnetic term).It should be noted that small-scale ionospheric irregularities are not considered in this study due to the difficulty to model and characterize them by a deterministic model, although their effect may be significant in some space weather situations.

Ray tracing method
The ray tracing technique is commonly used for calculating the propagation path of an electromagnetic signal in a medium specified by a position-dependent refractive index, such as the Earth's atmosphere.It has become a significant tool for investigating GPS signals' propagation; particularly it has been used in GNSS RO technology to investigate how the ionosphere affects the accuracy of neutral atmospheric retrievals.It has also been used to validate how much the separation of the GPS dual-frequency signals contributes to the bending angle RIE (Syndergaard, 2000).Hoque andJakowski (2010, 2011) used this method to study the effects of the ionospheric bending and high-order ionospheric error terms on GNSS RO signals.Mannucci et al. (2011) used this method to analyze the magnitude of the bending angle RIEs under severe ionospheric storms.In this study, a 3-D numerical ray tracing technique was used to simulate GPS signals received by low earth orbit (LEO) satellites after propagating through the atmosphere-ionosphere system to realistically obtain bending angle profiles under various ionospheric conditions.

Simulation of realistic observations
Realistic excess phase observations can be simulated by superimposing RO measurement errors (e.g., receiver thermal noise, precise orbit determination (POD) error, local multipath and clock instability) onto the excess phase profiles produced from the ray tracing stage described above.According to the GRAS receiving system's error budgets and characteristics, we set the error parameters similar to those adopted in Steiner and Kirchengast (2005) and Liu et al. (2013).The receiver noise was modeled as thermal noise with a 150 K LEO antenna noise and a 10 Hz loop band width.The POD error was modeled as a kinematic POD error with along-ray velocity and acceleration errors of 0.05 mm s −1 and 0.05 µm s −2 , respectively.The radial position errors of the GPS and the LEO satellites were set to 0.2 and 0.4 m, respectively.The local multipath effect was modeled using a sinusoidal-shaped function with the multipath phase error amplitude and period set to 0.5 mm and 100 s, respectively.The modeling of the clock errors was based on a random walk model and the ground-based single-differencing clock correction method with the relative stability of the ground clock being set to a 1 s Allan deviation of 1 × 10 −13 .

Computation of bending angle RIEs
The following three-step procedure was performed to compute the bending angle RIEs: (1) simulating no-ionosphere bending angle profiles using MSIS-90 only (i.e., no ionospheric model was used in the ray tracing) for obtaining a reference; (2) simulating ionosphere-corrected bending angle profiles using both aforementioned neutral atmospheric and ionospheric models and performing the dual-frequency linear combination of Eq. (3); and (3) obtaining the bending angle RIE profiles by differencing the ionosphere-corrected and the no-ionosphere bending angle profiles obtained from the above two steps.Figure 2 presents the detailed flow chart of the simulation process.

Ensemble simulation scheme
To characterize and quantify daily-global-mean and dailyzonal-mean bending angle RIEs, the RO events occurring on 15 July 2008 were simulated using the European Meteorological Operational (MetOp) satellite as an example LEO (Edwards and Pawlak, 2000) satellite with an onboard GPS receiver for atmospheric sounding, GRAS (Silvestrin et al., 2000).Both rising and setting occultation events were used in our investigation.Quasi-realistic atmospheric and ionospheric states were simulated using the MSIS-90 and 3-D NeUoG models.The total number of the simulated RO events on the day was found to be 723, of which 26 extremely noisy events were regarded as outliers (the upper threshold for acceptable bending angle RIEs was set to 7 µrad).As a result, 697 events were used in our study.
Figure 3 shows the global distribution of these 697 events; the background color represents the vertical total electron content (VTEC) as a function of latitude and longitude at 12:00 UT on 15 July 2008 calculated by 3-D NeUoG under the medium solar activity conditions (1 VTEC unit = 10 16 electrons m −2 ).According to the global distribution of the VTEC, the six geographical zones that include one global and five latitudinal bands/regions were chosen as the data bins.It should be noted that in the equatorial daytime (EDT) region, 84 RO events occurring between 09:00 and 21:00 local time (LT) were selected.Since MetOp is a sun-synchronous satellite and its equator crossing times were 09:30 LT (descending note) and 21:30 LT (ascending note), most of the events in the EDT region occurred in the time periods of 08:00-11:00 and 20:00-23:00 LT, while no event occurred in the time periods of 00:00-07:00 and 12:00-19:00 LT (Pirscher et al., 2007;Foelsche et al., 2009).Hence, the bending angle RIEs and their statistics represent the RIE characteristics in the periods of 09:00-11:00 and 20:00-21:00 LT rather than the entire daytime.This is relevant particularly at low latitudes and less relevant at high latitudes.

Simulation cases
To investigate "pure" bending angle RIEs, realistic bending angle errors and the effects of the ionospheric spherical symmetry assumption on the RIEs, we modeled the ionosphere under the following three scenarios (for the aforementioned RO simulation stage (2) in Sect.2.2.1): (1) without the ionosphere (wi), (2) spherical symmetry of the ionosphere (ss) and (3) non-spherical symmetry of the ionosphere (ns).We also modeled observation conditions (for stage 4 in Sect.2.2.1) for perfect observation (op) and realistic observation (or) scenarios, which refer to simulated RO observations with and without observing system errors, respectively.
Hereafter, the abbreviations opwi, opss, opns, orwi, orss and orns, which are the combination of the above three types of ionospheric conditions and two types of observation system conditions, will be used to denote six cases of simulated RO data sets.these six cases, which are also simulated for three different solar activity levels for the four cases, including the ionosphere.For example, opwi refers to the RO observations simulated without the effects of observation system errors (op) and without the ionosphere (wi), while opns refers to RO observations simulated without the effect of observation system errors (op) and with non-spherical symmetry (ns) settings used for the NeUoG ionosphere model.As a sensitivity check, we compared the opwi bending angles with the bending angles derived from the refractivity of the MSIS-90 model by the Abel transform under the assumption of spherical symmetry.Both results showed negligible differences < 0.1 µrad for the purpose (Liu et al., 2013).Thus the opwi bending angle data set was used as the reference data set to calculate the bending angle differences of all the other five cases relative to this data set.For example, the opns data set's bending angle RIEs were calculated by differencing its bending angles with the co-located opwi bending angles.

Calculation of bending angle RIE statistics
The statistics of the bending angle RIEs were calculated by standard algorithms such as those summarized in Scherllin- Pirscher et al. (2011a).For a single RO event, its absolute bending angle RIE profile was calculated by where the two bending angle terms at the right-hand side were obtained by interpolation and h i denotes the impact height levels (i.e., the impact parameter minus the local radius of curvature).The no-ionosphere and ionospherecorrected bending angles, which were simulated with a realistic sampling rate of 50 Hz, were interpolated to a standard vertical grid with 1600 impact height levels (50 m spacing between levels) in the 10-90 km range, because the number of sampling points in each of the reference profiles in this height range was mostly in the range of 1500-1700.
To show relative bending angle RIE profiles, i.e., in the form of percentages, the following formula for relative bending angle RIEs r (h i ) was used to assess the effects of the  bending angle RIEs on high-altitude bending angle retrievals: Finally, the level-average bending angle RIE bias, standard deviation (SD, σ ) and 2σ confidence-level uncertainty (2σ un) of the level-average bending angle RIE bias for the ensemble of bending angle profiles in a region (such as the global region in the ensemble study) were calculated by where n is the total number of the RO events in the region and j denotes the individual events.

Results and discussion
Table 3 shows the parameters of location, direction, time and rising/setting flag of five representative RO events that were selected to display ionospheric cross-sectional views along the occulting ray paths.Occ.44 and Occ.648 are shown in Fig. 4 and the features of the bending angle RIEs of one exceptional RO event (Occ.530)compared against the other two more typical events (Occ.26 and Occ.631) are shown in Fig. 5.In Table 3, the azimuth column (relative to the northern direction, counter-clockwise) indicates the GPS-to-LEO ray path direction and the rising/setting column indicates the RO vertical scanning directions (i.e., from the Earth's surface to the atmospheric top or the other way around).From Fig. 4 one gets a helpful impression of the character of the asymmetries of the ionospheric electron density along the inbound and outbound segments of the ray paths and of the variation of vertical electron density with the increase of solar activity level.From (c) we see that, for all three events, the difference between α 1 and α 2 somewhat increases with increasing impact height.The maximum differences between α 1 and α 2 of the Occ.26 and Occ.613 events at the impact height of 80 km reach about 12 and 20 µrad, respectively.In addition, the wave-like curves of the α c profiles indicate that α c contains bending angle RIEs.Comparing Occ.530 with the two normal events, one can find that there are extremely large fluctuations in the Occ.530 event's α 1 , α 2 and α c profiles, meaning that the exceptional event contains very large bending angle RIEs; thus, this event should be excluded in the calculation of the ensemble statistics.Similarly, one can find the behaviors of the excess phases of these events in panel (a).The L1 and L2 excess phases of the three events are around −13 and −23 m, respectively, and after the ionospheric correction the L c excess phases are around 0 m, as should be the case.Panel (d) shows that the maximum and minimum bending angle RIEs of the exceptional event reach 19 and −19 µrad, respectively, and the standard deviation of the exceptional event's bending angle RIEs is 5.7 µrad, which is about 20 times that of the normal events.This magnitude of bending angle RIEs is exceptional and so is regarded as an outlier.In this study, before performing statistical calculation, a value of 7 µrad was used as the threshold of outliers.Similarly, the maximum excess phase RIE of larger than 200 mm and the standard deviation of the exceptional event's excess phase RIEs of about 10 times the standard deviations of the normal events were found.We note that the approach in Liu et al. (2013) was used to calculate the standard deviations of the excess phases and bending angle RIEs for Fig. 5.The causes of the 26 outlier profiles rejected from the ensemble in total are investigated in a separate study; they may partly be physical (anomalous asymmetry effects) and partly technical (small discontinuities in the NeUoG refractivity field model perturbing the ray tracing in rare cases).

Ensemble simulation results
For our analyses of ensemble bending angle RIEs, the aforementioned six simulation cases combined with the six data zones were used to generate 36 ensemble bending angle data sets.Except for the opwi and orwi cases, each of the remaining 24 data sets includes results simulated under the three solar activity levels (F10.7 = 70, 140, 210).Of the 36 data sets, the six data sets GLO opwi, NHH opwi, NHM opwi, EDT opwi, SHM opwi and SHH opwi were used as the reference bending angle data sets for their corresponding zones in order to calculate the ensemble bending angle biases and their statistics according to Sect.2.3.3.
Figures 6-9 show the ionosphere-corrected bending angle profiles and their statistical RIE results (biases, standard deviations, 2σ confidence level uncertainties) for several representative data sets.The mean reference bending angle profiles in the left column refer to the mean of no-ionosphere bending angles that is from the corresponding opwi data set.In all these figures, the UM (upper mesosphere), LM (lower mesosphere), US (upper stratosphere) and LS (lower stratosphere) refer to the four impact height layers of 65-80, 50-65, 35-50 and 15-35 km, respectively.
Figure 10 illustrates relative bending angle RIEs and their statistics for the GLO opss and GLO orns representative data sets.The absolute and relative bending angle RIE biases were calculated by Eqs. ( 5) and ( 6), respectively, and their statistics by Eqs. ( 7)-( 9).The RIE bias profiles were calculated by averaging the ensemble bending angle RIE biases.The bias's 2σ confidence level uncertainties were calculated by Eq. ( 9), which is also effectively a 95 % confidence level of the ensemble mean bending angle RIEs.
Tables 4 to 7 present the layer-average absolute and relative bending angle RIE biases, standard deviations and their 2σ confidence level uncertainties in the following four impact height layers: the UM, LM, US and 20-35 km (which is a partial range of the LS) for the 24 data sets simulated with the ionospheric effects included.In these tables the layeraverage bending angle RIE biases and standard deviations of the four atmospheric layers were calculated by averaging the level-average bending angle RIE biases and their standard deviations over the 300 individual impact height levels in each impact height layer (50 m level spacing over 15 km).In fact, the layer-average bias is the average of biases of all the levels in the layer, and the layer-average standard deviation is approximately equal to the standard deviation of all the level RIE values in the same layer.
Since the correlation of bending angle errors at a given impact height level with their lower and upper neighboring levels are small, and the error correlation height range is approximately 1 km (Steiner and Kirchengast, 2005), selecting one effective sampling point in every 1 km height range should be sufficient for the calculation of a layer's 2σ confidence  level uncertainty using Eq. ( 9).In other words, about 15 independent effective sampling points should be used to calculate the layer's 2σ confidence level uncertainty.For example, when using Eq. ( 9) to calculate a layer's 2σ confidence level uncertainty, for a layer height extent of 15 km and the total number of RO events in the data set being N, the total number of independent effective sampling points in the layer is n = N × 15.Tables 4 to 7 list the corresponding values of n for each of the geographic zones.This use of 15 effective sampling points instead of the 300 individual (but highly correlated) levels is essential in order to not overestimate the smoothing of noise that is possible based on the data set.
Comparing the statistical results among the four impact height layers in each of Tables 4 to 7, we can see that, generally, the higher the impact height layer, the larger the layeraverage absolute bending angle RIE bias and the smaller the standard deviation.However, the relative bias results show that both layer-average biases and standard deviations increase with increasing impact height.Comparing among the three columns in each table for the same data set but different solar activity levels, one can see that, at the same impact height layer, the layer-average bending angle RIE bias and standard deviation increase with increasing solar activity level.The same results can also be seen from the comparison among the three panels in each column in Figs. 6 to 10.In Tables 4 to 7, the layer-average bending angle RIE biases are in the magnitude of less than 1 µrad.In terms of sign, the small RIE biases have a clear tendency to be negative, which is in line with the results of other studies in which real RO observation data were used (Sokolovskiy et al., 2009;Danzer et al., 2013).

Bending angle RIEs without observing system errors
To investigate the magnitude and characteristics of the pure bending angle RIEs (i.e., without observing system errors), a perfect observing system was used to simulate the opss and opns global ensemble bending angles.These results are shown in Fig. 6 and Tables 4-5.As shown in Fig. 6, the global-mean reference bending angle decreases exponentially with impact height, with values of about 1660 µrad at the impact height of 20 km, about 140 µrad at 35 km, about 16 µrad near the stratopause (about 50 km) and about 0.1 µrad at the mesopause (80 km).The global-mean bending angle RIE biases of both GLO opss and GLO opns data sets fluctuate around 0 and exhibit an obvious negative tendency.We can also see that the standard deviations of the level-mean bending angle RIE biases of the GLO opss data set are smaller than those of GLO opns due to the assumption of ionospheric spherical symmetry.In fact, the ionosphere is not spherical symmetric in reality, and hence the results of the GLO opns data set are closer to actual global bending angle RIEs of real data.
From Table 5 one can see that the layer-average absolute bending angle RIE biases of the GLO opns data set in the US at low, medium and high solar activity levels are −0.004,−0.007 and −0.015 µrad, respectively, and their corresponding standard deviations are 0.33, 0.41 and 0.53 µrad.In the LM, their mean bias values are −0.006,−0.008 and −0.02 µrad, and their corresponding standard deviations are 0.30, 0.39 and 0.52 µrad.The corresponding layer-average relative bending angle RIE biases in the US are 0.00, −0.02 and −0.03 % and in the LM are 0.01, −0.18 and −0.37 %.Evidently this reflects a clear increase with increasing solar activity.
Comparing layer-average bending angle RIE biases and their standard deviations in the same impact height layer at the same solar activity level but from the two different data sets, GLO opss (Table 4) and GLO opns (Table 5), we can see the effects of the ionospheric spherical symmetry assumption on the RIEs.The layer-average bending angle RIE biases in either US or LM impact height layer and at each of the three activity levels from the two data sets are very close.This suggests that the effects of the ionospheric spherical symmetry assumption on the layer-average bending angle RIE biases are fairly small.However, a comparison of the standard deviations of the layer-average absolute bending angle RIE biases at the same solar activity level and in the same impact height layer between the GLO opss and GLO opns data sets illus-trates that the assumption of the ionospheric spherical symmetry results in clear decreases in the standard deviations at the low, medium and high solar activity levels in the US by 18, 24 and 17 %, respectively, and in the LM by 20, 28 and 25 %, respectively.
According to Figs. 3 and 4, the EDT zone has the maximum VTEC and the steepest gradient of ionospheric electron density.As a result, the RO events in the EDT data zone have larger bending angle RIE biases than the other four zones, whilst their standard deviations at high solar activity level are smaller than those of the SHM data zone.As shown in Table 5 and Fig. 7, the layer-average bending angle RIE biases of the 84 events in the EDT opns data set in the UM, LM, US and the 20-35 km impact height layers under low solar activity level are −0.014,−0.009, −0.006 and −0.011 µrad, respectively.Their standard deviations are 0.30, 0.34, 0.37 events occurred in the LT periods of either 09:00-11:00 or 20:00-21:00, during which the ionospheric electron density values are small; otherwise, the bending angle RIE biases and their standard deviations would be even larger.Therefore, the results in this study reflect the characterizations of the bending angle RIEs of the sun-synchronous MetOp/GRAS LEO satellite system; fully local-time covering systems would also experience sometimes higher RIEs depending on local-time conditions.
From Table 5 and Fig. 7 we can see that the SHM opns data set contains 137 RO events.The layer-average bending angle RIE biases in the UM, LM and US under the low solar activity level are −0.002,−0.006 and −0.004 µrad, respectively, and their corresponding standard deviations are 0.32, 0.34 and 0.35 µrad; under the medium solar activity level, the layer-average bending angle RIE biases values are −0.015,−0.010 and −0.011 µrad, respectively, and their corresponding standard deviations are 0.38, 0.38 and 0.42 µrad; under the high solar activity, the mean bending angle RIE biases are −0.034,−0.023 and −0.024 µrad, respectively, and their corresponding standard deviations are 0.73, 0.78 and 0.76 µrad.This implies that the SHM zone experiences higher noise lev- Comparing the statistics of the bending angle RIEs of the EDT opns and SHM opns data sets in Table 5, we can find that the layer-average bending angle RIE biases of SHM are smaller.This is due to the reduction in its VTEC and electron density gradient.In contrast, the standard deviations of the layer-average RIE biases of the SHM opns data set are larger, especially under the high solar activity level, due to the following two factors.First, the RO events occurred in the whole data zone, including daytime and nighttime events, since these were all selected in the SHM opns data set.Additionally, there were some events in the LT period of 11:00-20:00 in the middle-latitude data set; thus the differences among the bending angle RIEs of its RO events are larger than those of the EDT data zone.Secondly, in the middle data zone, the inbound and outbound ray paths of the latitudinal direction RO events experience ionospheric asymmetry when passing through the ionosphere and the bending angle RIEs of these RO events are larger, whereas the longitudinal direction RO events' RIEs are smaller; therefore the difference between the latitudinal and longitudinal direction RO events' bending angle RIEs are larger than those in the EDT data zone.
Comparing the statistics of the bending angle RIEs at the same solar activity level in the same impact height layer but from the SHM opns and NHM opns data sets in Table 5, we can see that the layer-average RIE biases of the NHM data set are larger due to the equatorial anomaly in the North-ern Hemisphere.The standard deviations under the low and medium solar activity levels from both data sets are almost the same, but under the high solar activity level the RIEs of the NHM are smaller.
Similarly, from Table 5, a comparison of the SHH opns and NHH opns data sets shows that the layer-average absolute bending angle RIE biases and their standard deviations in SHH are smaller.In contrast, the layer-average relative bending angle RIE biases of SHH opns are larger due to the fact that the magnitude of the zonal-mean bending angle in SHH is smaller.Obviously, in all the six data sets, SHH exhibits the smallest layer-average absolute bending angle RIE biases and standard deviations.Furthermore, the differences of the layer-average bending angle RIE biases and their standard deviations resulting from the three solar activity levels are very small.This implies that the winter hemispheric highlatitudinal zone does not suffer more from high solar activity.The reason is that most of the area in that zone is perpetual night and the ionosphere is less strongly ionized.However, to add a remark of caution, the high-latitude errors are expected to generally be somewhat underestimated in the present endto-end simulations since the NeUoG model does not account for auroral zone and polar cap ionization anomalies that especially under high solar activity levels would render the highlatitude ionosphere much less smooth than represented by the NeUoG model.

Realistic bending angle RIEs including observing system errors
In order to first investigate the influence of the pure observing system errors (OSEs) on the global-mean and zonal-mean bending angle profiles, we simulated six orwi data sets.The resulting errors and their statistics were also calculated for the GLO, EDT and SHM data zones using Eqs.( 5)-( 9).The results in these three zones are shown in Fig. 8 and indicate that the bending angle residual biases from OSEs in all four atmospheric layers are essentially 0 and that the associated standard deviations in the US and LM layers amount to around 0.4 µrad.
As defined in Sect.2.3.2, the orns simulation cases, which were simulated by superimposing various observing system errors, can reflect the magnitude and characteristics of realistic bending angle RIEs (including residual observing system errors) of the MetOp/GRAS RO observing system.As shown in Fig. 9 and Table 7, the layer-average absolute bending angle RIE biases in the US at the low, medium and high solar activity levels are −0.003,−0.008 and −0.015 µrad, respectively.Their corresponding standard deviations are 0.51, 0.56 and 0.66 µrad.In the LM, they are −0.004,−0.012 and −0.021 µrad, and their corresponding standard deviations are 0.49, 0.55 and 0.64 µrad.From the same data set, the layeraverage relative bending angle RIE biases under the three solar activity levels in the US are −0.01,−0.02 and −0.03 %, and in the LM they are −0.01,−0.22 and −0.37 %.The con-  A comparison of the GLO orss (Fig. 9 and Table 6) and GLO orns (Fig. 9 and Table 7) data sets shows the effects of the ionospheric spherical symmetry assumption on the bending angle RIEs under the realistic observation system scenarios.The assumption of ionospheric spherical symmetry leads to reductions in the absolute RIE standard deviations under low, medium and high solar activity levels in the US by 6, 11 and 11 %, and in the LM by 6, 13 and 14 %, respectively.
Comparing the bending angle RIE biases and their standard deviations in the same impact height layer at the same solar activity level but from two data sets GLO opns (Fig. 6 and Table 5) and GLO orns (Fig. 9 and Table 7), one can see that both the bending angle RIE and the residual error from the OSEs are relevant error sources in the US, LM and UM, but with different characteristics.For example, in terms of layer-average bending angle RIE biases, the values from GLO opns are about equal to those of GLO orns in each of the three layers at each of the three solar activity levels.This suggests that the negative mean biases in the global-mean bending angle profiles in the US, LM and the UM are mainly caused by the bending angle RIEs, since the residual OSEs are random errors that produce zero biases.This can be confirmed by the fact that in the SHH orns and orss data sets all the layer-average bending angle RIE biases in the US, LM and UM at all the three solar activity levels are around 0. The bending angle RIEs in these two data sets are determined by the OSEs due to the fact that most area in the winter hemispheric high-latitudinal zone was under low ionization.
In terms of the mean standard deviations, in the US for example, the "orns-opns" ratios of the global ensemble bending angle RIEs under low, medium and high solar activity levels are about 1.54, 1.32 and 1.25, respectively; in the LM their values are 1.64, 1.43 and 1.23; and in the UM their values are 1.65, 1.43 and 1.25.This leads to the conclusion that the residual observing system errors can make the standard deviations of the total residual errors increase by about 23 to 64 % so that OSEs are clearly a main contributor to the residual random errors.A comparison between the GLO orss and GLO opss data sets under the assumption of ionospheric spherical symmetry also shows similar results.

Overall discussion
To provide a convenient summary view of the effects of the bending angle RIEs and the residual errors from OSEs on realistic bending angle errors, we selected the UM layer to show this; the other three layers (not shown) exhibit similar RIE characteristics.The UM layer-average bias and standard deviation of the bending angle RIEs for each of the 697 RO events in the GLO f70opns, GLO f210opns, GLO orwi, GLO f70orns and GLO f210orns data sets were calculated using the approach in Liu et al. (2013).Figure 11 shows a histogram distribution of the number of RO events as a function of the UM layer's bias (left column) and standard deviation (right column) of the bending angle errors for each of the five data sets.In the panels, the red lines denote the values of layer-average biases and standard deviations from Tables 5 and 7; the blue lines denote the event-average biases and standard deviations calculated by averaging the layeraverage biases and standard deviations of each of the 697 events in the data set; the green lines denote the median biases and the median standard deviations of all the 697 RO events.
As shown in Fig. 11, in terms of the bending angle bias, all the 697 RO events in the GLO f70opns, GLO opwi and GLO f70orns data sets spread in the ranges of −0.13 to 0.09, −0.11 to 0.13 and −0.13 to 0.13 µrad, respectively.In both GLO f210opns and GLO f210orns, some of the bias values are beyond the range of −0.15 to 0.15 µrad due to the high solar activity level.For all the five data sets, most of the RO events fall into the range −0.03 to 0.03 µrad, and the −0.03 to 0.01 bins have the largest number of RO events.Obviously, the histogram of the GLO opwi data set is symmetric with the layer-average, event-average and median bending angle RIE bias of 0. This implies that the residual OSEs have a Gaussian distribution.The distributions of the other four histograms are skewed to the left so they have negative layer-average (equal to event-average) and negative median biases.The layer-average bending angle biases of GLO f70opns, GLO f210opns, GLO f70orns and GLO f210orns are −0.004,−0.019, −0.003 and −0.019 µrad, respectively.These further confirm that the negative biases in the dailyglobal-mean bending angle profiles are mainly caused by the bending angle RIEs.

Summary and conclusions
Previous studies have proven the feasibility of using the GNSS RO technique for observing decadal and longer period climate in the UTLS height region.However, errors in atmospheric RO observations due to ionospheric influences could be considerable, and thus the ionospheric error effects on atmospheric variable retrievals need to be mitigated for highaccuracy operational weather forecasting and climate monitoring.In addition, according to the GCOS (2010) report, upper-air temperature is one of the essential climate variables to be observed systematically for climate monitoring.However, RO-retrieved temperature in the US and higher layers cannot yet meet the accuracy required.Due to the bending and separation of GPS dual-frequency ray paths when passing through the atmosphere, RIEs in bending angles are one of the main error sources of bending angles in the MS and US layers even after the standard dual-frequency ionospheric correction of bending angles is performed.These bending angle RIEs will propagate into refractivity, temperature and other related RO retrievals.It is therefore essential to characterize and quantify the bending angle RIEs for improved error mitigation algorithms, which was the focus of this study.
We have quantified the global-mean and zonal-mean bending angle RIE biases in the US, LM and UM height layers using realistic end-to-end simulations.The characteristics and quantities of the ensemble bending angle RIEs in the three layers at three different solar activity levels in various observational conditions and latitudinal zones have been investigated and analyzed using six RO simulation cases (opwi, opss, opns, orwi, orss, orns; Table 2) studied in six geographic data zones (GLO, NHH, NHM, EDT, SHM, SHH; Table 1).
The results show that the mean bending angle RIE biases in the US, LM and UM are small but have a negative tendency.In climatologies of retrieved profiles, these biases may lead to systematic errors.The magnitude of the bending angle RIE biases and their standard deviations increase with in- creasing solar activity level.With regard to the regional characteristics and quantities of the bending angle RIE statistics, the EDT data zone shows the largest RIE biases.The layeraverage bending angle RIE biases in the two middle-latitude data zones are slight smaller than those in the EDT, while their standard deviations are slightly larger.
Generally high-latitude data zones have the smallest bending angle RIE biases and standard deviations, which is partly also an effect of disregarding in the large-scale ionospheric modeling used the medium-and small-scale ionization variability present in the auroral zone and polar cap ionosphere.The differences in the bending angle RIE biases and their standard deviations between the SHH and NHH data zones are noticeable.The layer-average bending angle RIE biases in SHH in all four height layers under all three solar activity levels and in all RO simulation cases are close to 0. This is due to the weakly ionized ionosphere during the polar night in the winter hemisphere.The layer-average bending angle biases in the NHH data zone (polar summer) are also small but still noticeable, especially under the middle and high solar activity levels.
Compared to RIE studies using real RO data, such as the recent one by Danzer et al. (2013), the magnitude of the neg-ative biases found from our simulations is roughly only half the one estimated from real data.One reason may be that our simulations are based on a large-scale ionospheric model, which is somewhat realistic but still does not represent all effects of the real ionosphere (e.g., possible contributions to RIE biases from asymmetric small-scale irregularities or from the geomagnetic higher-order refraction term are not included).However, the analyses of real data may in part attribute bias contributions from non-RIE observational errors into the RIE estimations (e.g., from residual clock correction or multipath errors).Either way, further work to reconcile current bias estimations from simulations and real-data analyses will be needed.
Given the estimated magnitude of the bending angle RIE biases in the US, LM and UM in tropical and mid-latitude areas, it is clearly desirable to develop algorithms for mitigating the RIEs for high-accuracy weather forecasting and climate monitoring applications.Fortunately, residual bending angle errors induced by the observing system errors can be considered essentially random errors, so that the biases can be interpreted as largely of ionospheric origin, while the bending angle standard deviations contain a main contribution also of observation system errors.Therefore, the bending angle RIE biases can possibly be partially mitigated under realistic observation conditions, if the determining factors of the RIE, such as solar activity, latitudinal zone and time are known.Alternatively, a correction may be possible on a perevent basis, using knowledge of the ionospheric conditions in the ionospheric inbound and outbound regions of the events.
Overall, the characterization and quantification of the bending angle RIEs performed in this study contributed essential and valuable insights into RIE characteristics to aid modeling of RIEs and improve atmospheric profile retrievals in the future.The results help in this sense to inform future RIE mitigation schemes that will improve upon the use of the standard linear ionospheric correction of bending angles and will also provide explicit uncertainty estimates, all for the benefit of further improved quality and climate utility of RO data over the stratosphere and mesosphere.

Figure 3 .
Figure 3. Distribution of the 697 RO events used in the statistical analysis: the upward-pointing triangles denote rising events, and the downward-pointing triangles denote setting events; the five latitudinal zones (NHH, NHM, SHM, SHH, EDT; Table 1) are divided by latitude-circle lines (solid lines); in the EDT region, although all the RO events during the day and night are shown here, only those RO events occurring during the local daytime (LT) 09:00-21:00 belong to the region.The background color map illustrates the ionospheric VTEC (1 VTEC unit = 10 16 electrons m −2 ) at 12:00 UTC on 15 July 2008 under the medium solar activity level (F10.7 = 140).

Figure 4
Figure4depicts the vertical electron density (1 electron density unit (EDU) = 10 11 electrons m −3 ) distribution along the Occ.44 and Occ.648 occultation event planes (i.e., latitudinal direction versus altitude plane at the time 11:00 LT) at three solar activity levels (F10.7 = 70, 140, 210).Three representative ray paths are shown for each event.From Fig.4one gets a helpful impression of the character of the asymmetries of the ionospheric electron density along the inbound and outbound segments of the ray paths and of the variation of vertical electron density with the increase of solar activity level.

Figure
Figure 5a-d show a comparison of excess phases (in meters), excess phase RIEs (in unit of mm), bending angles and bending angle RIEs (both in unit of µrad), respectively, in the impact height range of 40-80 km at medium solar activity level.An exceptionally noisy RO event (Occ.530) and two typical events (Occ.26 and Occ.631) are illustrated.Panels (a) and (c) present the three events' L1, L2 and ionosphere-corrected excess phases and bending angles.From (c) we see that, for all three events, the difference between α 1 and α 2 somewhat increases with increasing impact height.The maximum differences between α 1 and α 2 of the Occ.26 and Occ.613 events at the impact height of 80 km reach about 12 and 20 µrad, respectively.In addition, the wave-like curves of the α c profiles indicate that α c contains bending angle RIEs.Comparing Occ.530 with the two normal events, one can find that there are extremely large fluctuations in the Occ.530 event's α 1 , α 2 and α c profiles, meaning that the exceptional event contains very large bending angle RIEs; thus, this event should be excluded in the calculation of the ensemble statistics.Similarly, one can find the behaviors of the excess phases of these events in panel (a).The L1 and L2 excess phases of the three events are around −13 and

Figure 5 .
Figure 5.Comparison of excess phases (a), excess phase RIEs (b), bending angles (c) and bending angle RIEs (d), respectively, in the impact height range of 40-80 km at the medium solar activity level between the exceptional RO event Occ.530 (red) and two typical RO events, Occ.26 (blue) and Occ.631 (green), respectively.

Figure 6 .
Figure6.Ionosphere-corrected bending angle profiles and their statistical results for the global ensemble "opss" data set (left) and "opns" data set (right), respectively.In both composite panels, the left column illustrates the ionosphere-corrected bending angle profiles (mean reference bending angle (blue line), ensemble bending angles (grey lines) and their mean bending angle (MBA) in the UM (orange), LM (red), US (cyan) and LS (green)) under low (top), medium (middle) and high (bottom) solar activity levels.The right column shows their corresponding statistical RIE results including the bias (thin lines), standard deviation (thick lines) and the bias' 95 % confidence level uncertainty (dark green lines).

Figure 7 .
Figure 7. Ionosphere-corrected bending angle profiles and their statistical results for the EDT (left) and SHM (right) "opns" data sets.The figure layout and legends are the same as in Fig. 6.

Figure 8 .
Figure 8. Ionosphere-corrected bending angle profiles and their statistical results for the "orwi" data set.The figure layout and legends are the same as in Fig. 6.

Figure 9 .
Figure 9. Ionosphere-corrected bending angle profiles and their statistical results for the GLO "orss" (left) and "orns" (right) data sets.The figure layout and legends are the same as in Fig. 6.

Figure 10 .
Figure10.Relative RIEs and their statistics for the GLO ensemble "opss" (left column) and "orns" (right column) data sets.The bias (thin lines), standard deviation (thick lines) and the bias' 95 % confidence level uncertainty (dark green lines) are depicted.

Figure 11 .
Figure11.Relationship between the histogram distribution of the number of RO events and the bias (left panels) and standard deviation (right panels) of the bending angle RIEs in the UM from the GLO f70opns, GLO f210opns, GLO opwi, GLO f70orns and GLO f210orns data sets (top to bottom).The red, blue and green lines denote layer average, event average and median biases (left panels) and standard deviations (right panels) of bending angle RIEs, respectively.

Liu et al.: Quantifying residual ionospheric errors in bending angles
Flow chart of the RO end-to-end simulation process for bending angle RIEs; for a description see Sect.2.2.

Table 1 .
Table 2 summarizes the detailed definition of Geographic zone definitions.

Table 2 .
Definition of end-to-end simulation cases for various ionospheric conditions and solar activity levels under the assumptions of no observing system errors or realistic observing system errors.

Table 3 .
Parameters of the RO example events illustrated in Figs. 4 and 5.

Table 5 .
Absolute and relative biases, including their 2σ (95 %) confidence ranges, and absolute and relative standard deviations of bending angle RIEs in the same layout as in Table4but for the opns case.

Table 6 .
Absolute and relative biases, including their 2σ (95 %) confidence ranges, and absolute and relative standard deviations of bending angle RIEs in the same layout as in Table4but for the orss case.Bending angle residual ionospheric error (RIE) estimates for realistic observing system and spherical symmetry case (orss)

Table 7 .
Absolute and relative biases, including their 2σ (95 %) confidence ranges, and absolute and relative standard deviations of bending angle RIEs in the same layout as in Table4but for the orns case.Bending angle residual ionospheric error (RIE) estimates for realistic observing system and nonspherical symmetry case (orns)