Improved model for correcting the ionospheric impact on bending angle in radio occultation measurements

The standard approach to remove the effects of the ionosphere from neutral atmosphere GPS radio occultation measurements is to estimate a corrected bending angle from a combination of the L1 and L2 bending angles. This approach is known to result in systematic errors and an extension has been proposed to the standard ionospheric correction that is dependent on the squared L1/L2 bending angle difference and a scaling term (κ). The variation of κ with height, time, season, location 10 and solar activity (i.e. the f10.7 flux) has been investigated by applying a 1D bending angle operator to electron density profiles provided by a monthly median ionospheric climatology model. As expected, the residual bending angle is well correlated (negatively) with the vertical TEC. However, κ is more strongly dependent on the solar zenith angle. Furthermore, over the height range of interest (40-80 km) κ is approximately linear with height. A simple κ model has also been developed. It is independent of ionospheric measurements, but incorporates geophysical dependencies (i.e. solar zenith angle, solar flux, 15 altitude). The global mean error (i.e. bias) and the standard deviation of the residual errors are reduced to -2.2×10 rad and 2.0×10 rad respectively. Although a fixed scalar κ also reduces bias for the global average the selected value of κ (14 rad) is only appropriate for a small band of locations around the solar terminator. In the day time, the scalar κ is consistently too high and this results in an over correction of the bending angles and a positive bending angle bias. Similarly, in the night time, the scalar κ is too low. However, in this case, the bending angles are already small and the impact of the choice of κ is less 20 pronounced.


Introduction
It has been demonstrated that, by using variational data assimilation techniques, GPS radio occultation (GPS-RO) measurements can be assimilated into operational numerical weather prediction (NWP) systems to improve the accuracy of temperatures in the upper troposphere and lower-middle stratosphere (Healy and Thépaut, 2006;Poli et al., 2009;Rennie, 2010). In particular, GPS-RO measurements reduce stratospheric temperature biases in NWP systems and this indicates that such measurements could have an increasingly important role in climate monitoring and climate reanalyses (Poli et al., 2010;Steiner et al., 2013). Notwithstanding the benefits of GPS-RO for the neutral atmosphere, it remains necessary to consider the effect of the ionosphere on the measurements.
Vorob 'ev and Krasil'nikova (1994) (hereafter referred to as VK94) proposed a method of combining the GPS-RO bending angles measured at two frequencies (L1 and L2) to provide a first-order correction for the ionosphere. VK94 also showed that the first-order correction leaves a systematic bending angle bias that increases as a function of the electron density squared, integrated over the vertical profile. The relationship between the bias and electron density suggests that the bending angle biases should vary diurnally and as a func- tion of the 11-year solar cycle. This has been demonstrated by various authors, e.g. Kursinski et al. (1997), Mannucci et al. (2011) and Danzer et al. (2013). Healy and Culverwell (2015) have proposed a modification to the standard bending angle correction to reduce the residual systematic ionospheric errors. The modification introduces a new second-order term that is a function of the square of L1 and L2 bending angle difference and a weighting term (κ). The aim of this work is to investigate the variation of κ with height, time, season, location and solar activity (i.e. the F10.7 flux). This has been done by applying a 1-D bending angle operator to electron density profiles provided by the NeQuick monthly median ionospheric climatology model (Nava et al., 2008). As well as examining the variations in κ, a κ model has been developed. It is independent of ionospheric measurements and therefore simple to implement in an operational system, but does incorporate the relevant geophysical dependencies (i.e. solar zenith angle, solar flux). The expectation is that, since NeQuick is a reasonable median model of the ionosphere, the κ model derived from it will also exhibit reasonable statistics, though this has not been proven.
Radio occultations, the VK94 ionospheric correction procedure and the proposed modified correction are described in Sect. 2. Examples of how κ varies with height, location and solar activity are presented in Sect. 3. Models for κ are proposed and assessed in Sect. 4 and the discussion and conclusions are given in Sects. 5 and 6.
2 Radio occultation and ionospheric corrections Hardy et al. (1994), Kursinski et al. (1997) and Hajj et al. (2002) provide a comprehensive description of the GPS-RO technique. In summary, the GPS satellites transmit on two L-band channels (L1, L2) at f 1 = 1575.42 MHz and f 2 = 1227.60 MHz and the signals are received by a satellite in low earth orbit (LEO) (Fig. 1).
The standard approach (Abel transform) for inverting GPS-RO measurements requires the assumption of spherical symmetry. Under that assumption, the bending angle of the ray between the GPS satellite and a receiver in LEO is where i = 1, 2 depending on the frequency, a is the impact parameter, r t is the tangent height of the ray path and n i is the refractive index. The impact parameter is given by Horizontal gradients will result in residual errors in the inversion. However, it is expected that these errors are random; therefore, they should not affect monthly or seasonal climatologies.
To a first-order approximation, the refractive index comprises terms dependent on the neutral atmosphere refractivity (N n ), the ionospheric electron density (n e ) and the frequency (f ) squared: where k = 40.3 m 3 s −2 . Therefore, the measured L1 and L2 bending angles are different from each other and contain both neutral and ionospheric components. The standard approach taken in operational RO processing centres is to estimate a corrected neutral atmosphere bending angle (α c ) using the VK94 approach: where the L1 and L2 bending angles (α L1 and α L2 respectively) are interpolated to a common impact parameter. The first-order approximation neglects terms involving higher powers of the frequency and the earth's magnetic field; however, these have little effect on the residual bending angle errors (Syndergaard, 2000). One benefit of this approach is that it is based on the standard parameters estimated by the RO retrieval system and does not require a priori information about the ionosphere. One downside is that a systematic bending angle error remains (see Eq. 5 of Healy and Culverwell, 2015). The bending angle error has a dependence on the electron density squared, which indicates that it will vary with the solar cycle. This has been recognised as a potential source of bias in climatology products (Danzer et al., 2013). Healy and Culverwell (2015) have proposed a modification to the standard ionospheric correction of the following form: Day of month Not used The day of month is used to linearly interpolate between two monthly coefficient files. This prevents step changes in electron density at month boundaries.
hmE Hard coded to 120 km Hard coded to 110 km. This is a more reasonable value. However, a more sophisticated model will be implemented in future; i.e. Chu et al., 2009. Bottom-side taper Displays a discontinuity at 90 km that can produce artefacts in bending angle estimations.
Bottom-side taper added using a tanh function.
where the κ term compensates for the systematic residual error in the standard approach. An appropriate value for κ has been investigated using simple analytic functions for the ionosphere  and using a ray tracer through a 3-D ionospheric model (Danzer et al., 2015), though it should be noted that this study was limited to a low latitude band because of noise in the simulation system. It has been shown that κ generally falls in the range of 10 to 20 rad −1 and a simple scalar model, κ ∼ 14, provides a good first approximation, improving the accuracy of the "dry" temperature retrievals (Danzer et al., 2015). Nevertheless, it is clear that κ will vary as a function of height, local time, season, location and solar activity and therefore it is possible that existing ionospheric climatology models could be used to compute an improved correction term by modelling the monthly mean, temporal and spatial variations of κ more realistically.

Examples of κ dependencies
A monthly median 3-D ionospheric model (in this case NeQuick) and a 1-D bending angle operator (based on Eq. 1) can be used to estimate the residual ionospheric error and thereby estimate values for κ.

NeQuick
NeQuick is a monthly median ionospheric electron density model developed at the Aeronomy and Radiopropaga-  (Nava et al., 2008). The model is based on the Di Giovanni-Radicella (DGR) model (Di Giovanni and Radicella, 1990), which was modified for the European Cooperation in Science and Technology (COST) Action 238 to provide electron densities from ground to 1000 km. The model has been designed to have continuously integrable vertical profiles which allows for rapid calculation of the total electron content (TEC) for transionospheric propagation applications. The current version of NeQuick can be run up to a height of 20 000 km and a  variant is used in the Galileo Global Navigation Satellite System (GNSS) to calculate ionospheric corrections (Angrisano et al., 2013). NeQuick is a "profiler" which makes use of three profile anchor points at the E layer peak, the F1 peak and the F2 peak. To specify the anchor points it uses the layer critical frequencies (foE, foF1, foF2) and the F2 maximum usable frequency factor (M3000(F2)) (Davies, 1965). foE is determined using a solar zenith angle model, foF1 is assumed to be proportional to foE during daytime and zero during nighttime, and foF2 and M3000(F2) are derived from the ITU-R (CCIR) coefficients in the same way as the International Reference Ionosphere (IRI) (Bilitza and Reinisch 2008).
Between 100 km and the peak of the F2 layer, NeQuick uses an electron density profile based on the superposition of five semi-Epstein layers (Epstein, 1930;Rawer, 1983); i.e. the Epstein layers have different thickness parameters for their top and bottom sides. The top side of NeQuick is a simplified approximation to a diffusive equilibrium. A semi-Epstein layer represents the model top side with a heightdependent thickness parameter that has been empirically determined.
The model used in this work is the University of Birmingham's translation of the NeQuick v2.0.2 from FORTRAN into Python. Very minor (negligible) differences in results are observed due to the use of different interpolation routines. The Python code has been largely vectorised to increase the speed of operation. Some additional modifications have been made and are described in Table 1.

κ estimation
In each of the examples shown in the following sections the same basic procedure has been followed to estimate the value of κ: 1. Use NeQuick to estimate a vertical profile of electron density.
2. Convert the electron density (n e ) to the refractive index (n i ) using the first-order approximation n i = 1 − 40.3 n e /f 2 i for each frequency (L1 and L2). 3. Estimate bending angle using the 1-D observation operator for L1 and L2.
Atmos. Meas. Tech., 11, 2213-2224, 2018 www.atmos-meas-tech.net/11/2213/2018/  Since no neutral atmosphere is included in the estimate of the refractive index, α c should be zero if VK94 provides a perfect correction. Any non-zero values are representative of the residual ionospheric error ( α), which, from Eq. (5), is modelled as Since the bending angles are known, this can be rearranged to provide an estimate of κ as a function of the impact parameter: In real data the corrected bending angles increase rapidly towards the surface. This means that the impact of any residual error becomes less insignificant below approximately 40 km. Furthermore, the VK94 correction assumes that the ray impact parameter/tangent height is below the ionosphere (i.e. the electron density is zero). Consequently, the main area of interest for κ estimation is between 40 and 80 km. It is in this region where the residual error from the ionospheric correction is likely to be a major contributor to the overall error budget of neutral atmosphere retrievals.

Height dependence
The Figs. 2 to 5 show two examples of the vertical electron density profile, the L1 / L2 bending angles, the residual error and κ. The test parameters are given in Table 2. Over the height range of interest (40-80 km), Fig. 5 shows that κ is approximately linear with tangent height, but its gradient is dependent on the local time.

Geographic dependence
The geographic dependence of bending angle correction can be demonstrated by plotting maps of the TEC (Fig. 6), residual bending angle (Fig. 7) and κ (Fig. 8). In this case, the test parameters are given in Table 3. As expected, the residual bending angle is well correlated (negatively) with the verti-    cal TEC. However, κ is more strongly dependent on the solar zenith angle, indicating that the TEC-dependent component of the residual error is largely modelled by the squared L1 / L2 bending angle difference term in the correction and that κ is modelling other features such as changes in hmF2.

Solar cycle dependence
The solar cycle dependence of κ has been investigated by estimating κ at a tangent height of 60 km above London for each day over the last 60 years (Table 4). The results (Fig. 9) show that κ is negatively correlated with F10.7; i.e. κ is low when the vertical TEC is large, which occurs when F10.7 is high. Furthermore, the dynamic range of κ is considerably smaller than that of the F10.7 (and hence TEC and bending angle), varying by a factor of approximately 50 % compared to approximately 300 % for F10.7. This, again, is indicative of the TEC-dependent component of the residual error being largely modelled by the squared L1 / L2 bending angle difference term in the correction.
4 Models of κ

Introduction
Section 3 has presented examples of how κ can vary spatially and with solar cycle. In this section, simple models of κ will be assessed in order to evaluate their potential to reduce the residual bending angle errors in the VK94 correction. Three models will be considered:  κ equals zero (zero-κ): this represents the current situation with the unmodified VK94 correction.
κ is a scalar (scal-κ): this is the approach proposed by Healy and Culverwell (2015).
κ is a function of latitude, longitude, solar zenith angle and solar flux (func-κ).
In order to build the models a set of 25 000 κ estimates were generated from NeQuick using random drivers (uniformly distributed over the ranges in Table 5). The true solar flux is used for each randomly selected day and year. A further independent set of 25 000 κ estimates was also generated using the same random parameter ranges to act as a test data set.

Scalar κ
The random κ values are shown in Fig. 10. The median value is marked by the horizontal line and has value of 14 rad −1 . This value is used as the scalar model.
accounted for by the scalar model. Figure 8 is suggestive that κ is a function of solar zenith angle -this is a convenient parameter to use since it embodies the position, local time and season. Figures 11, 12 and 13 show κ as a function of solar zenith angle, F10.7 and altitude respectively. Note that the solar zenith angle has been extended to π radians to account for when the sun is below the horizon. The figures indicate broadly linear dependencies in all cases; therefore, the following model is proposed: where f 10.7 is the F10.7 flux (sfu, 1 sfu = 10 −22 Wm −2 Hz −1 ), χ is the solar zenith angle (rad) and h is the height above the ground (km); a, b, c and e are scalars to be found by fitting the model to the data.
The Python code curve_fit from the scipy.optimize package has been used to fit the model. The parameter results and the associated variances are shown in Table 6. A plot of the NeQuick estimated κ compared to the func-κ is shown in Fig. 14. Figure 15 shows the geographic distribution of funcκ at 12:00 UT in June and December at 60 km altitude and with an F10.7 of 150. These maps can be directly compared with those in Fig. 8.

Bending angle error reduction
The second set of 25 000 randomly distributed points has been used to assess the reduction in residual bending angle for each of the κ models (zero-κ, scal-κ and func-κ). Figure 16 shows a histogram of the residual bending angle errors for the full data set. The bending angle error statistics are in Table 7.
Both the scal-κ and func-κ results are an improvement over the zero-κ results. In the case of the scal-κ, both the  standard deviation and the mean error (i.e. bias) of the residual bending angle errors are reduced by an order of magnitude compared to the zero-κ results (from −1.3 × 10 −8 and 2.2 × 10 −8 for the zero-κ case to 5.4 × 10 −9 rad and 1.5 × 10 −9 rad respectively; Table 7). In the case of the func-κ, the standard deviation and the mean error (i.e. bias) of the residual bending angle errors are further reduced to 2.0 × 10 −9 rad and −2.2 × 10 −10 rad respectively. Although the scal-κ reduces bias for the global average, the geographic variation of κ (shown in Figs. 8 and 15) makes it clear that the selected value of κ (14 rad −1 ) is, in fact, only appropriate for a small band of locations around the solar terminator. The effect of this is clear if the residual error statistics are considered for daytime and nighttime separately. Figures 17 and 18 show histograms for residual bending angle for day and night respectively. In the daytime, the scalκ is consistently too high and this results in an overcorrection of the bending angles and a positive bending angle bias (Table 7). Similarly, in the nighttime, scal-κ is too low and this results in a negative bending angle bias. However, in this  case, the bending angles are already small and the impact of the choice of κ is less pronounced (Table 7). Conversely, the func-κ results in a negative bending angle bias in the daytime and positive bias at night. In both cases, the bending angle biases are significantly lower than those produced with scal-κ.

Discussion
Many studies of ionospheric refraction of transionospheric radio waves have shown that, in addition to the level of ionisation, the shape of the vertical electron density profile plays a significant role, e.g. Jakowski et al. (1994) and Jakowski (2008, 2010). It is important to remember that the functional model of κ has been created by fitting κ derived from NeQuick. NeQuick is based on the standard CCIR databases of foF2, foE and M3000F2 and therefore provides a reasonable median model of the F and E regions' peaks; however, it is not certain that NeQuick is a good median representation of the layer shapes. Furthermore, the κ model is derived from 1-D estimates of the bending angle and so does not take nonspherical structures into consideration. The approach, therefore, has been to model kappa with minimal complexity to avoid a close fitting to NeQuick that may be inappropriate in reality. Additional terms have been also trialled in the model (such as local time), but these have not shown any significant improvement of the model presented in this paper. Given the limitation of the ionospheric model and the bending angle estimation, the results are indicative that a simple kappa model may be used, but further testing with real data must be done to validate this.

Conclusions
Using the random selection of vertical profiles from the NeQuick the median κ has been shown to be 14 rad −1 and this is therefore an appropriate value for κ if it is to be represented by a single scalar. This value agrees well with the result from Healy and Culverwell (2015) and is in the range suggested by Danzer et al. (2015). Representing κ as a scalar has the advantage of simplicity and is appropriate if climate reprocessing centres are focused on ensuring that global average biases are removed. However, it has been demonstrated that such an approach can lead to significant differences in the residual bending angle bias between day and night. In the day, the results indicate that the bending angle bias switches sign from −3.3 × 10 −8 rad for no correction to +7.6 × 10 −9 rad for the scalar κ correction.
This limitation can be overcome using the simple κ function model. This approach does not require independent ionospheric measurements and so remains easy to implement. It should be noted that the κ model is based on a monthly median ionospheric model. Whilst this is a starting point it will be necessary to work with climate reprocessing centres to develop an effective validation strategies of the bending angle corrections. It would also be useful to assess the sensitivity of stratospheric climatologies to the bending angle bias and standard deviation bounds determined by this study. Furthermore, the magnitude of other error terms (i.e. non-symmetry; Zeng et al., 2016) should be assessed in light of these results.
Data availability. The NeQuick variant used to produce the results in this paper can be requested from the Abdus Salam International Centre for Theoretical Physics (ICTP), Trieste, Italy.