A modification to the standard ionospheric correction method used in GPS radio occultation

A modification to the standard bending-angle correction used in GPS radio occultation (GPS-RO) is proposed. The modified approach should reduce systematic residual ionospheric errors in GPS radio occultation climatologies. A new second-order term is introduced in order to account for a known source of systematic error, which is generally neglected. The new term has the form κ(a)× (αL1(a)− αL2(a)) , where a is the impact parameter and (αL1, αL2) are the L1 and L2 bending angles, respectively. The variable κ is a weak function of the impact parameter, a, but it does depend on a priori ionospheric information. The theoretical basis of the new term is examined. The sensitivity of κ to the assumed ionospheric parameters is investigated in one-dimensional simulations, and it is shown that κ ' 10– 20 rad. We note that the current implicit assumption is κ = 0, and this is probably adequate for numerical weather prediction applications. However, the uncertainty in κ should be included in the uncertainty estimates for the geophysical climatologies produced from GPS-RO measurements. The limitations of the new ionospheric correction when applied to CHAMP (Challenging Minisatellite Payload) measurements are noted. These arise because of the assumption that the refractive index is unity at the satellite, made when deriving bending angles from the Doppler shift values.


Introduction
GPS radio occultation (GPS-RO) measurements are now routinely assimilated into operational numerical weather prediction (NWP) systems (Healy and Thépaut, 2006;Cucurull et al., 2007;Aparicio and Deblonde, 2008;Poli et al., 2009;Rennie, 2010).Using variational data assimilation techniques, it has been shown that they provide accurate temperature information in the upper troposphere and lower/middle stratosphere, complementing the information provided by satellite radiance measurements.Specifically, it has been shown that GPS-RO measurements reduce stratospheric temperature biases in NWP systems.The value of GPS-RO in climate reanalyses has been demonstrated (Poli et al., 2010), and a number of studies outlining climate-monitoring applications have been reported (e.g.Leroy et al., 2006).A particularly noteworthy activity is the "RoTrends project", where the major GPS-RO processing centres have processed common data sets in order to estimate structural uncertainty in the geophysical retrievals (Ho et al., 2009(Ho et al., , 2012;;Steiner et al., 2013).
Overall, the published results indicate that GPS-RO could have an increasingly important role in climate monitoring in the coming years, particularly in the stratosphere.However, one area of potential concern that could affect climate-monitoring applications is the impact of "residual ionospheric errors" on the geophysical retrievals.These arise because the measurement is sensitive to both the neutral atmosphere and the ionosphere.The first-order "ionospheric correction" commonly used in the GPS-RO processing (Vorob'ev and Krasil'nikova, 1994) produces systematic retrieval errors that will vary as a function of the 11-year solar cycle (Danzer et al., 2013).
The impact of residual ionospheric errors on GPS-RO retrieval accuracy has been discussed by a number of authors.Vorob'ev and Krasil'nikova (1994) were the first to propose correcting for the ionosphere at the bending-angle level rather than phase delays.They derived an integral expres-Published by Copernicus Publications on behalf of the European Geosciences Union.
sion to estimate a source of systematic residual error in their bending-angle approach, showing that these errors increased as a function of the electron density squared, integrated over the vertical profile.They also argued that horizontal gradients only weakly affected the accuracy of their method, as the first-order impact of the gradients is removed in the correction procedure.Kursinski et al. (1997, Table 2) estimated the accuracy of the bending-angle method as a function of solar and diurnal cycles, based on simulations with onedimensional (1-D), spherically symmetric Chapman layer ionospheres.The largest error was at daytime solar maximum conditions, reaching −6.5 % (around −0.3 µrad) at 60 km.Syndergaard (2000) introduced an improved phase correction method with similar error characteristics to the bendingangle method.A "major" and "minor" error term was identified, and it was noted that the bending-angle correction removed the major term, but the removal of the minor term required additional a priori information.Interestingly, Syndergaard argued that in practice the term in the ionospheric refractive index related to the earth's magnetic field could be neglected when assessing the residual ionospheric errors in GPS-RO applications.More recently, Mannucci et al. (2011) investigated the errors in simulations based on threedimensional ionospheric fields from an ionospheric assimilation model.Based on these results they suggested that ionospheric errors at a height of 20 km are too large for climatemonitoring applications, during daytime at solar maximum conditions.Furthermore, they emphasised the need for improved ionospheric correction schemes to mitigate this problem.Danzer et al. (2013) investigated an improved ionospheric correction based on the statistical differences of measurements taken at night and day in order to reduce systematic geophysical retrieval errors.Complementary to Danzer et al. (2013), Liu et al. (2013Liu et al. ( , 2015) ) performed end-to-end simulations to investigate residual ionospheric errors remaining after the dual-frequency bending-angle correction.They also found clear evidence for negative residual bending-angle bias, albeit smaller in size than the observation-based studies.
In this study, we present a new, relatively simple approach for reducing the systematic residual ionospheric error originally identified and investigated by Vorob'ev and Krasil'nikova (1994).The method still requires a priori ionospheric information, which is combined with the observed bending-angle information in order to estimate a correction term.This work has possible relevance to the generation of accurate monthly mean geophysical climatologies of the stratosphere, retrieved from GPS-RO measurements.In particular, it may be useful when producing GPS-RO climatologies with the new average-bending-angle method, which has been discussed recently (Ao et al., 2012;Gleisner and Healy, 2013;Danzer et al., 2014).The new method does not correct the residual errors arising from horizontal gradients in the ionosphere, or those caused by the earth's geomagnetic field.We expect that these effects are likely to produce random noise on individual profiles, but to first order they should av-erage out, and not affect monthly and seasonal climatologies.However, this has not been proven here.
The theoretical basis of the new method for reducing the systematic residual ionospheric errors will be outlined, and we will demonstrate the method in simulations with one-dimensional, spherically symmetric model ionospheres, where the electron number density is only a function of the vertical co-ordinate.More detailed simulations which expand on these results -based on modelling through complex, three-dimensional ionospheres -will be reported by Danzer et al. (2015).
In Sect.2, the ionospheric correction procedure proposed by Vorob'ev and Krasil'nikova (1994), which has been adopted at most processing centres, will be described and the new model for the residual error will be introduced.Onedimensional simulations testing the new approach are presented in Sect.3. The discussion and conclusions are given in Sect.4, where the strengths and limitations of the new approach are outlined.

Modified ionospheric correction of GPS-RO bending angles
It has been demonstrated that GPS-RO measurements provide useful temperature information in the stratosphere.However, GPS-RO is not a direct measurement of temperature profile information, and a retrieval system is required to estimate the geophysical variables such as pressure, temperature and geopotential height.The basic components of the GPS-RO geophysical retrieval system are outlined in Kursinski et al. (1997) and Hajj et al. (2002).We will focus on the ionospheric correction step here.
The GPS satellites transmit signals at two L band frequencies, f 1 = 1575.42MHz and f 2 = 1227.60MHz.The removal of the ionospheric bending to first order is possible because the ionosphere is dispersive at the GPS L band frequencies.Assuming spherical symmetry, the bending angle, α Li , at impact parameter, a, is where i = (1, 2), depending on the frequency; r t is the tangent height of the ray path; and n i is the refractive index.The refractive index is approximated with where N n is the neutral refractivity, the constant k 4 = 40.3m 3 s −2 , and n e (r) is the electron number density.This expression neglects higher-order terms involving higher powers of the frequency f i and the earth's magnetic field, but Syndergaard (2000) has found that this has no appreciable impact on the residual bending-angle errors.
As noted earlier, the ionospheric correction is usually performed at the bending-angle level at most operational processing centres, based on the approach suggested by Vorobev and Krasilnikova (1994) (VK94, hereafter).The "corrected" neutral atmosphere bending angle, α c , at impact parameter a, is estimated with where α L1 and α L2 are the bending angles for L1 and L2 signals, respectively, interpolated to a common impact parameter value, a.A common impact parameter is equivalent to a common ray tangent point height when the ray-path tangent point is below the ionosphere.The residual errors generally increase as the ionospheric electron densities increase (e.g.Kursinski et al., 1997;Syndergaard, 2000;Mannucci et al., 2011;Danzer et al., 2013;Liu et al., 2013).
One of the strengths of Eq. ( 3) is that it is based on simple physics, and it does not require a priori information about the ionospheric state.It is generally accepted that ionospheric correction schemes that improve upon this will require some a priori ionospheric state information (e.g.Syndergaard, 2000).
It is interesting to note that VK94 provided an estimate of a systematic residual ionospheric bending-angle error, although -with the exception of Syndergaard (2000) and Danzer et al. (2013) -this does not seem to have received much attention (see Eq. 22, VK94).Neglecting the neutral contribution to the refractive index in order to simplify the mathematics (n i = 1−k 4 n e (r)/f 2 i ), VK94 asssumes that the ray impact parameter/tangent height is below the ionosphere, so that n e (r t ) = 0 and a = r t .They insert n −1 i 1 + k 4 n e (r)/f 2 i and n 2 i 1 − 2k 4 n e (r)/f 2 i into Eq. 1 and expand the denominator of the integral to obtain an approximation that is second order in n e , Applying the standard ionospheric correction (Eq. 3) to the second-order approximation (Eq.4) produces a systematic residual term given by 1 (r 2 − a 2 ) 3/2 dr.
(5) 1 The numerator of VK94 (Eq.22) differs slightly.It has a factor (3r 2 −2a 2 ) in Eq. ( 22), whereas we find (2r 2 −a 2 ).This difference has no significant impact on the magnitude of the error estimate, but Fig. 1 suggests our expression is correct.This error term arises even for the simplest case of a spherically symmetric ionosphere, with no magnetic field.Integration by parts shows that the integral in Eq. ( 5) is always positive (see Eq. A4), so this error biases the corrected bending angles negative, meaning the values produced by Eq. ( 3) are consistently too small.The error occurs because the L1 and L2 signals have different ray paths.The dependence on n 2 e indicates that the bias will depend on the solar cycle.The present study is concerned with estimating and correcting this specific source of bias in the bending angles, in order to improve GPS-RO geophysical climatologies.However, we emphasise that this requires some a priori assumptions about the ionospheric state.
In this study, we propose a modification to the standard ionospheric correction of the form where the new term, κ(a)×(α L1 (a)−α L2 (a)) 2 , compensates for the systematic residual error identified in Eq. ( 5).The physical justification for the new term and estimates of κ(a) based on both simulations through one-dimensional, spherically symmetric ionospheres, and analytical solutions of the bending-angle equation are presented in Sect.3.

Results
We have confirmed the accuracy of Eq. ( 5) for onedimensional, spherically symmetric ionospheres, where the electron number density is only a function of the vertical coordinate and neutral atmosphere refractivity index is unity.Figure 1 shows an example where L1 and L2 bending angles have been computed by integrating Eq. ( 1) for Chapman layer ionosphere.The Chapman layer electron number density profile is given by The Chapman layer peaks at r m = R e + 300 km, where R e is the radius of the earth.It has a width of H = 75 km and a peak electron number density of n max e = 3 × 10 12 m −3 .These parameters are most relevant to daytime, solar maximum conditions.The computed bending angles are then corrected with Eq. ( 3), but since there is no neutral bending in the simulation the departure from zero is the computed ionospheric residual.The analytic estimate of the residual error using Eq. ( 5) is also shown.It is clear that the agreement between computed and estimated errors below ∼ 70 km is extremely good.The divergence above 80 km is because of the assumption in the analytical expression that the electron number density is zero at the ray tangent point.The residual errors are actually small when compared with the L1 and L2 bending angles.For example, at 60 km α L1 = 215 µrad and α L2 = 354 µrad.The residual error at 60 km is 0.27 µrad, of order 0.1 % of α L1 .However, the neutral bending angle falls exponentially with height, and at 60 km it has a climatological average value of 4-5 so the bias in the corrected value is of order ∼ 6 % (Kursinski et al., 1997).Note that the magnitude of the residual error is comparable to estimates presented in Mannucci et al. (2011), and by Danzer et al. (2013) at solar maximum conditions.Consequently, we believe Eq. ( 5) is a significant component of the total bias in the corrected bending angles.
There is a linear relationship between the residual error value, α(a), and (α L1 (a) − α L2 (a)) 2 when we adjust the Chapman layer peak electron density, n max e , keeping both the height, r m , and width, H , values constant.Inspection of the ionospheric bending-angle integral (Eq. 1) and error term (Eq.5) reveals they scale as n max e and (n max e ) 2 , respectively.This relationship suggests that the sensitivity of the systematic residual error with respect to the peak electron density, n max e , can be modelled in terms of (α for the Chapman layer used in Fig. 1. κ is 15.8 rad −1 at the surface, falling to 11.6 rad −1 at 100 km.The corrected bending angles increase exponentially towards the surface, meaning the impact of the residual error becomes less significant there, so we are most interested in the κ values around ∼ 40 km and above.In Appendix A, we provide analytical expressions for κ(a), based on three 1-D model electron density profiles where both the ionospheric bending angles and residual error term can be approximated analytically.This approach provides a theoretical basis for Fig. 2 and gives insight into how the κ(a) scales with the assumed ionospheric parameters.The model ionospheres, shown in Fig. 3, are (1) an idealised Chapman layer; (2) a slab profile centred at r m ; and (3) an asymmetric triangle, peaking at r m .These analytical forms and the computed value (Eq.8) are shown in Fig. 4. The model ionospheres have the same vertically integrated total electron content (τ e ) value and n max e values.The agreement between the computed κ (Eq.8) and the analytical approximations is encouraging, with all values in the range of ∼ 10-20 rad −1 near 50 km.The analytical approximation for the idealised Chapman layer (Eq.A13) gives κ values that are ∼ 20 % too low compared to the computed Chapman values.This is because the bending-angle and error integrals are approximated and simplified in order to provide analytical solutions (Eqs.A7 and A8).However, Eq. (A13) reproduces the relatively weak dependence on a, and it suggests that the sensitivity to the assumed r m and H is approximately κ ∝ (r m −a) 1/2 /H .The slab model κ values are higher than the other models, but this is related to the atypical "shape factor" of the slab ionosphere, which is defined as The shape factors for a Chapman layer, asymmetric triangle and slab models are 0.66, 0.67 and 1.0, respectively.As a consequence, the slab model κ will not tend to the same limit as the Chapman layer for small H , unless an adjusted shape factor is arbitrarily introduced.The asymmetric triangle gives the best agreement with the computed values; this is because the integrals (Eqs.A3 and A4) can be evaluated accurately, and the shape factor is reasonable.
The sensitivity of computed κ(a) values to the assumed ionospheric parameters is illustrated in Figs. 5 and 6.The width of the Chapman layer is 50, 75 and 100 km.The height of the peak electron number density is 250, 300 and 400 km.These results are reasonably consistent with the sensitivity suggested by the analytical expressions, showing that κ falls as H is increased and increases when r m is increased.The shape of κ differs when H = 100 km (Fig. 5) and r m = R e +250 km (Fig. 6) because of the non-zero electron density values at the ray tangent height.

Discussion and conclusions
This study has focused on a source of systematic residual ionospheric error, originally noted in VK94.The results suggest it may be beneficial to introduce an additional term in the standard bending-angle ionospheric correction scheme (Eq.6).
This introduces a new term, κ(a)×(α L1 (a)−α L2 (a)) 2 .The variable κ is a weak function of the impact parameter, a, but it does depend on an ionospheric state model through an assumed width parameter, H , and the height of the peak electron density, r m .We have used both simulations and analytical approximations to explore the sensitivity of κ(a) to these parameters, and the range of these results suggests κ 10-20 rad −1 is a reasonable approximation.
We note that the modified ionospheric correction scheme will be important in the generation of GPS-RO geophysical climatologies if the temporal variability of κ(a) within the 11-year solar cycle is small compared to the variability of the (α L1 (a) − α L2 (a)) 2 term.This hypothesis has been investigated recently by Danzer et al. (2015).They have performed simulations with a three-dimensional ionospheric model, for a period spanning an 11-year solar cycle.They  8), and the analytical estimates for an idealised Chapman layer (Eq.A13), a slab ionosphere (Eq.A14) and an asymmetric triangle (Eq.A17).The profiles have the same n max e and vertical TEC values.have found reasonable agreement between the theoretical κ values and those derived from the simulations, and that the use of suitably chosen values of κ can significantly reduce errors in upper-stratospheric climatological temperatures derived from GPS-RO data.Nevertheless, further testing with other three-dimensional ionospheric models should be undertaken.
It must be accepted that κ(a) still depends on an underlying ionospheric model and that it may be difficult to validate the new correction scheme directly against independent observations.However, we have replaced the sensitivity of the residual error with respect to the peak electron density with a sensitivity with respect to the measured L1 and L2 bending-angle values.Nevertheless, there will still be some uncertainty in κ, and this will project into the uncertainty in the geophysical climatologies.We note that the current bending-angle ionospheric correction scheme (Eq. 3) implicitly assumes κ = 0.This is probably adequate for NWP applications, but neglecting the second-order term in this way will introduce time-varying biases in GPS-RO climatologies.
A realistic, non-zero κ(a) of the correct order of magnitude should reduce these biases, even if it does not remove them completely.Furthermore, it may be possible to obtain more accurate κ(a) values from ionospheric climatology models, or reanalyses based on ionospheric data assimilation systems.This should be investigated in future work.In any case, the sensitivity of the GPS-RO geophysical climatologies to the assumed κ(a) will be of interest to users, and it should be a component of the climatology uncertainty estimation.This uncertainty is not currently captured in the RoTrends project (Steiner et al., 2013) because all the processing centres currently assume κ = 0.
One limitation of the new ionospheric correction approach that should be highlighted is regarding the application to CHAMP (Challenging Minisatellite Payload) data.The bending angles we receive are not direct measurements but are derived from a Doppler shift value assuming the refractive index at the low-earth-orbit (LEO) satellite is unity, meaning the electron density at the LEO is assumed to be zero.This assumption introduces systematic errors in both the L1 and L2 bending-angle values which scale as n e (r leo )/f 2 i (Hajj and Romans, 1998), where r leo is the position of the LEO satellite.These errors are large for CHAMP measurements because the altitude is ∼ 420 km, but they do not affect the corrected bending angles because they cancel out in the standard linear ionospheric correction scheme (Eq.3).Unfortunately, these errors do not cancel when computing (α L1 − α L2 ) 2 , and they can introduce an error of ∼ 60 % in this term for CHAMP measurements at solar maximum conditions.The error is about an order of magnitude smaller for COSMIC and GRAS measurements (∼ 5 %) because the altitude is ∼ 800 km.However, clearly this limitation should be considered when combining CHAMP with other measurements into a single time series.
In summary, we have investigated a systematic residual error in the standard GPS-RO ionospheric correction.A modification to the standard ionospheric correction has been suggested which may be particularly important when generating geophysical climatologies from GPS-RO measurements.The initial results are promising, and they suggest the modified approach should be considered at the GPS-RO processing centres.

Appendix A: Analytical approximations for κ(a)
Simple 1-D ionosphere models are useful for understanding the magnitude of κ(a) and the sensitivity to the assumed ionospheric parameters, H and r m .The approach is to approximate the linearised ionospheric contribution to the bending integral, (i = 1, 2) and the systematic residual error term, assuming that the Chapman layer is sufficiently peaked at r m .This assumption appears to be acceptable because it is made when approximating both the bending and error integrals, and the final result depends on the ratio of these expressions.The solutions are where H is the Chapman layer width (Eq.7) and τ e = n max e H √ 2π e is the vertical total electron content.We can substitute Eq. (A7) into Eq.(A8): We then assume the ionospheric bending, α i , can be approximated with the measurements (α L1 , α L2 ), composed of both ionospheric and neutral bending, and the ionospheric corrected value, α c (a), giving and The residual error is then (A12) Finally using Eq. ( 8), Note the scaling with 1/H and r m , which define the assumed ionospheric model.

A2 Slab ionosphere
A similar procedure can be followed for a slab ionosphere, where the electron density is a constant value over a vertical interval of total width 2H , centred on r m , and is zero elsewhere.The bending angle and error integrals (Eqs.A3 and

Figure 1 .
Figure 1.Comparing analytically estimated residual ionospheric errors (Eq.5) and computed values, for a spherically symmetric Chapman layer ionosphere, with no neutral bending.The Chapman layer peaks at 300 km, has a width of H = 75 km and a peak electron number density of n max e = 3 × 10 12 m −3 .The computed errors are found by integrating Eq. (1) for both L1 and L2, and then correcting the L1 and L2 bending angles with Eq. (3).

Figure 2 .
Figure 2. The computed κ (Eq.8) for a Chapman layer which peaks at 300 km has a width of H = 75 km and a peak electron number density of n max e = 3 × 10 12 m −3 .

Figure 3 .
Figure 3.The profiles of the electron number density used to provide analytical estimates of κ.

Figure 4 .
Figure 4.The computed κ given by Eq. (8), and the analytical estimates for an idealised Chapman layer (Eq.A13), a slab ionosphere (Eq.A14) and an asymmetric triangle (Eq.A17).The profiles have the same n max

Figure 5 .
Figure 5.The impact of the assumed width, H = (100, 75, 50) km, on the computed κ values for a Chapman layer.

Figure 6 .
Figure 6.The sensitivity of the computed κ to the assumed r m .The altitudes, z m = r m −R e , of the peak electron density tested are 250, 300, and 400 km.
)analytically.We then substitute the ionospheric bendingangle solution into the error equation, and then rearrange to find an equation of the form−κ(a) × (α L1 − α L2 ) 2 .If the electron density is zero at the tangent height, n e (a) = 0, Eqs.(A1) and (A2) can be integrated by parts to remove the electron density gradient terms, giving ) and (A4) cannot be integrated analytically for a Chapman layer.Hence, the ionospheric bending-angle and residual error integrals are approximated with