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

Introduction Conclusions References Tables Figures


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 Figures Back Close Full 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 datasets in order to estimate structural uncertainty in the geophysical retrievals (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 phasedelays.They derived an integral expression 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 one-dimensional, 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 bending angle 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 ar-Figures

Back Close
Full gued 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 three-dimensional 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 climate monitoring 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.
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 climatogies with the new average bending angle method, 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 average 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).Introduction

Conclusions References
Tables Figures

Back Close
Full 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.One dimensional 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.

Ionospheric correction of GPS-RO measurement
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., Kursinki et al., 1997;Syndergaard, 2000;Mannucci et al., 2011;Danzer 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 to simplify the mathematics, their correction method (Eq. 3) assumes the ionospheric bending is well ap-Introduction

Conclusions References
Tables Figures

Back Close
Full This assumes that the L1 and L2 bending does not lead to significant deviations from the straight line paths.VK94 set n e to 0 at r t , which therefore equals a, and expand the denominator of the left hand side of Eq. ( 4) to estimate the second-order term for both the L1 and L2 frequencies.Applying Eq. ( 3) to the expanded, 2nd order bending angle equations produces a systematic residual term given by1 , 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.Introduction

Conclusions References
Tables Figures

Back Close
Full

Results
We have confirmed the accuracy of Eq. ( 5) for one-dimensional (1-D), spherically symmetric ionospheres, where the electron number density is only a function of the vertical co-ordinate, 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 These parameters are most relevant to day time, 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 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 µrad, 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 (α L1 (a) − α L2 (a)) 2 .Figure 2 shows 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 The agreement between the computed κ (Eq.7) 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 1185 Introduction

Conclusions References
Tables Figures

Back Close
Full 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.68, 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, and this is because the integrals (Eqs.A3 and A4) can be evaluated acurately, 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: The variable κ is a weak function of 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 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 climatogies.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 obtained 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 data.The bending angles we receive are not direct measurements, but are derived from a Doppler shift value assuming the refractive index at the 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 Introduction

Conclusions References
Tables Figures

Back Close
Full approximate the linearised ionospheric contribution to the bending integral with, ) and the systematic residual error term, analytically.We then substitute the ionospheric bending angle solution into the error 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.6) and TEC = n and The residual error is then Finally using Eq. ( 7), Note the scaling with 1/H and r m , which define the assumed ionospheric model.

A2 Slab ionosphere 10
A similar procedure can be followed for a slab ionosphere, where the electron density is 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 A4) can be solved more Introduction

Conclusions References
Tables Figures

Back Close
Full easily in this case.The κ value is where l = (r m − a)/H, and, provided we approximate r 2 − a 2 by (r m + a)(r − a) in Eqs. ( A3) and (A4).

A3 Asymmetric triangle
The asymmetric triangle ionosphere has a peak electron density at r m .The electron density falls with height above and below the peak value, and is zero at (r m − H 1 ) and (r m + H 2 ), and is zero elsewhere.The ratio (H 2 /H 1 ) 2.152 is chosen so that the fraction of the TEC above r m is the same as with the Chapman layer ionosphere.Again Eqs.(A3) and (A4) can be evaluated relatively easily.Approximating r 2 − a 2 by (r m + a)(r − a) in Eqs.(A3) and (A4), we find, where l i = (r m − a)/H i .The functions are,

Conclusions References
Tables Figures

Back Close
Full  Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |proximated with an integral that is linear in n e , Discussion Paper | Discussion Paper | Discussion Paper | has a width of H = 75 km, and a peak electron number density of n Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 2 ,
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.7) are shown in Fig. 4. The model ionospheres have the same vertically integrated total electron content (TEC) value, and n max e values.The agreement between the computed κ (Eq.7) 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 Discussion Paper | Discussion Paper | Discussion Paper | 2 (α L1 (a) − α L2 (a)) + κ(a)(α L1 (a) − α L2 (a)) 2 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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, α i (a) 2a10 −6 k 4 and (A4) cannot be integrated analytically for a Chapman layer.Hence, the ionospheric bending angle and residual error integrals are approximated with, α i (a) 2a10 −6 k 4 f 2 then assume the ionospheric bending, α i , can be approximated with the measurements, (α L1 , α L2 ), composed of both ionospheric and neutral bending, and the iono-Discussion Paper | Discussion Paper | Discussion Paper | spheric corrected value, α c (a), giving Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | and B

Figure 1 .Figure 2 .Figure 3 .
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).