Annales Geophysicae Higher order ionospheric effects in GNSS positioning in the European region

After removal of the Selective Availability in 2000, the ionosphere became the dominant error source for Global Navigation Satellite Systems (GNSS), especially for the high-accuracy (cm-mm) demanding applications like the Precise Point Positioning (PPP) and Real Time Kinematic (RTK) positioning. The common practice of eliminating the ionospheric error, e.g. by the ionosphere free (IF) observable, which is a linear combination of observables on two frequencies such as GPS L1 and L2, accounts for about 99 % of the total ionospheric effect, known as the first order ionospheric effect (Ion1). The remaining 1 % residual range errors (RREs) in the IF observable are due to the higher – second and third, order ionospheric effects, Ion2 and Ion3, respectively. Both terms are related with the electron content along the signal path; moreover Ion2 term is associated with the influence of the geomagnetic field on the ionospheric refractive index and Ion3 with the ray bending effect of the ionosphere, which can cause significant deviation in the ray trajectory (due to strong electron density gradients in the ionosphere) such that the error contribution of Ion3 can exceed that of Ion2 (Kim and Tinin, 2007). The higher order error terms do not cancel out in the (first order) ionospherically corrected observable and as such, when not accounted for, they can degrade the accuracy of GNSS positioning, depending on the level of the solar activity and geomagnetic and ionospheric conditions (Hoque and Jakowski, 2007). Simulation results from early 1990s show that Ion2 and Ion3 would contribute to the ionospheric error budget by less than 1 % of the Ion1 term at GPS frequencies (Datta-Barua et al., 2008). Although the IF observable may provide sufficient accuracy for most GNSS applications, Correspondence to: Z. G. Elmas (isxzge1@nottingham.ac.uk) Ion2 and Ion3 need to be considered for higher accuracy demanding applications especially at times of higher solar activity. This paper investigates the higher order ionospheric effects (Ion2 and Ion3, however excluding the ray bending effects associated with Ion3) in the European region in the GNSS positioning considering the precise point positioning (PPP) method. For this purpose observations from four European tations were considered. These observations were taken in four time intervals corresponding to various geophysical conditions: the active and quiet periods of the solar cycle, 2001 and 2006, respectively, excluding the effects of disturbances in the geomagnetic field (i.e. geomagnetic storms), as well as the years of 2001 and 2003, this time including the impact of geomagnetic disturbances. The program RINEXHO (Marques et al., 2011) was used to calculate the magnitudes of Ion2 and Ion3 on the range measurements as well as the total electron content (TEC) observed on each receiver-satellite link. The program also corrects the GPS observation files for Ion2 and Ion3; thereafter it is possible to perform PPP with both the original and corrected GPS observation files to analyze the impact of the higher order ionospheric error terms excluding the ray bending effect which may become significant especially at low elevation angles (Ioannides and Strangeways, 2002) on the estimated station coordinates.

2000, the ionosphere became the dominant error source for Global Navigation Satellite Systems (GNSS), especially for the high-accuracy (cm-mm) demanding applications like the Precise Point Positioning (PPP) and Real Time Kinematic (RTK) positioning.
The common practice of eliminating the ionospheric error, e.g. by the ionosphere free (IF) observable, which is a linear combination of observables on two frequencies such as GPS L1 and L2, accounts for about 99 % of the total ionospheric effect, known as the first order ionospheric effect (Ion1).The remaining 1 % residual range errors (RREs) in the IF observable are due to the higher -second and third, order ionospheric effects, Ion2 and Ion3, respectively.Both terms are related with the electron content along the signal path; moreover Ion2 term is associated with the influence of the geomagnetic field on the ionospheric refractive index and Ion3 with the ray bending effect of the ionosphere, which can cause significant deviation in the ray trajectory (due to strong electron density gradients in the ionosphere) such that the error contribution of Ion3 can exceed that of Ion2 (Kim and Tinin, 2007).
The higher order error terms do not cancel out in the (first order) ionospherically corrected observable and as such, when not accounted for, they can degrade the accuracy of GNSS positioning, depending on the level of the solar activity and geomagnetic and ionospheric conditions (Hoque and Jakowski, 2007).Simulation results from early 1990s show that Ion2 and Ion3 would contribute to the ionospheric error budget by less than 1 % of the Ion1 term at GPS frequencies (Datta-Barua et al., 2008).Although the IF observable may provide sufficient accuracy for most GNSS applications, Correspondence to: Z. G. Elmas (isxzge1@nottingham.ac.uk)Ion2 and Ion3 need to be considered for higher accuracy demanding applications especially at times of higher solar activity.
This paper investigates the higher order ionospheric effects (Ion2 and Ion3, however excluding the ray bending effects associated with Ion3) in the European region in the GNSS positioning considering the precise point positioning (PPP) method.For this purpose observations from four European stations were considered.These observations were taken in four time intervals corresponding to various geophysical conditions: the active and quiet periods of the solar cycle, 2001 and 2006, respectively, excluding the effects of disturbances in the geomagnetic field (i.e.geomagnetic storms), as well as the years of 2001 and 2003, this time including the impact of geomagnetic disturbances.The program RINEX HO (Marques et al., 2011) was used to calculate the magnitudes of Ion2 and Ion3 on the range measurements as well as the total electron content (TEC) observed on each receiver-satellite link.The program also corrects the GPS observation files for Ion2 and Ion3; thereafter it is possible to perform PPP with both the original and corrected GPS observation files to analyze the impact of the higher order ionospheric error terms excluding the ray bending effect which may become significant especially at low elevation angles (Ioannides and Strangeways, 2002) on the estimated station coordinates.

Introduction
After removal of Selective Availability in 2000, the ionosphere became the dominant error source in Global Navigation Satellite Systems (GNSS) error budget (El-Rabbany, 2002).The ionosphere is a medium of free electrons and ions and as such its dispersive nature makes the ionospheric Published by Copernicus Publications on behalf of the European Geosciences Union.Z. G. Elmas et al.: Higher order ionospheric effects in GNSS positioning refractive index different from unity -this difference in the refractive index leads to the group and phase velocities of the GNSS signals to differ from each other during propagation through the ionosphere, such that the group velocity decreases (leading to the group delay i.e. code measurements longer than the geometric range) and the phase velocity increases (leading to the carrier phase advance i.e. carrier phase measurements shorter than the geometric range).
When the phase and group velocities are affected, the ray direction is also likely to be affected unless the wave is travelling perpendicularly to the gradient in ionospheric refractive index (Cairo and Cerisier, 1976).This effect, also known as the ray bending effect, is inversely proportional to the signal frequency and is highly dependent on the satellite elevation angle.The error due to the ray bending is orders of magnitude smaller than the first order ionospheric error; it is indeed comparable to that of the higher order error terms (Petrie et al., 2010).
The diffractive and ray-bending effects of the ionosphere on GNSS signals are neglected in this work; only the part of the Ion3 term quadratic in terms of the electron density is analysed, although the third order error due to ray bending may, under some conditions, exceed or be comparable in magnitude to the second order error term (Hoque and Jakowski, 2008).This negligence of the ray bending effect assumes that the GNSS signals travel along the straight LoS path between the receiver and satellite (Hoque and Jakowski, 2007) instead of two slightly different paths (bent and LoS paths).This assumption (of neglecting a bent path) then implies that the TEC and geomagnetic field effect along the signal path are the same for different (e.g.L1 and L2) signal frequencies.The RINEX HO program does not estimate corrections for the bending effect thus the corrected (for Ion3) measurements will still be contaminated by the ray bending effect.It should also be mentioned at this point that for the estimation of the Ion2 term, the ionosphere is assumed as a single thin shell at 450 km and the magnetic induction B is computed at the ionospheric pierce point (IPP) and not along the ray path of the signal.This assumption is expected to introduce some errors to the computation of the Ion2 term however it has less computational burden than the ray path approach.
Subsequently when positioning is performed with the corrected and uncorrected measurements an elevation cutoff angle of 10 • is considered in the PPP.As shown by Hoque and Jakowski (2007), the ray bending effect on the GNSS signals becomes significant especially at low satellite elevation angles -below 10 • ; thus this cutoff angle is thought to be an appropriate threshold for comparing the positioning results (considering the corrected and uncorrected observations) with negligible contribution by the ray bending effect.More detailed analyses about the impact of ray bending on the GNSS signal propagation and observations have been shown, among other researchers, by Hoque and Jakowski (2008), who provided an empirical formula for the geometric bending of the GPS signals; Hartmann and Leitinger (1984), who considered the geometric bending of the signals in their analysis of the RREs due to the atmosphere; and Petrie et al. (2010), who used the International Reference Ionosphere (IRI) 2007 model (Bilitza and Reinish, 2008) to estimate the potential size of the ray bending effect on the estimated GPS parameters and positioning.
The different orders of the ionospheric error terms (Ion1, Ion2 and Ion3, as denoted in this work for the first, second and third order terms) can have magnitudes that are observed to change according to the background solar, ionospheric and magnetic conditions.At times of high background solar activity, as during the peaks of the solar cycle or active days of an ionospheric storm, the greater amount of solar radiation causes increased levels of electron density in the ionosphere.This can cause the slant range delay on the GPS L1 signal link to be as large as 100 m in the uncorrected observable (the error contributed by the Ion1 term is about 10 to 100 m in general as shown by Klobuchar andKunches, 2003 andGrewal et al., 2007).In the ionospherically corrected (to the first order as in the dual frequency applications) observable, RREs can reach tens of centimetres (the Ion2 term contributes about several centimetres of range error as shown by Bassiri andHajj, 1992, 1993;Kedar et al., 2003;Fritsche et al., 2005;Hoque and Jakowski, 2006;Morton et al., 2009) and the Ion3 term about 1 cm or less, e.g. during disturbed ionospheric background conditions, which may happen due to geomagnetic storms, as discussed by Bassiri and Hajj (1993), Brunner and Gu (1991) and Kedar et al. (2003).
In general, most of the ionospheric range error can be eliminated depending on the method of positioning performed, i.e. stand-alone or differential.In stand-alone mode, users with a dual frequency receiver can account for the first order ionospheric effect by the IF observable, whereas users with a single frequency receiver can resort to an ionospheric model like the Klobuchar model (Leick, 2004).The IF observable can eliminate about 99 % of the total ionospheric effect (i.e. the Ion1 term) yielding an accuracy sufficient for most GNSS applications; an ionospheric model like the Klobuchar model, however, can correct about 50-60 % of the total ionospheric effect and gives limited performance for the users outside the mid-latitudes (Orus et al., 2002).In the differential mode, for short baselines, the ionospheric error can be eliminated by ionospheric corrections obtained from a reference station assuming a spatially and temporally correlated ionosphere (for such short baseline) between the user and reference.However, during adverse ionospheric conditions spatial and temporal correlation of the ionospheric errors can decrease.
For high accuracy demanding GNSS applications like PPP and RTK, especially during the peaks of the solar cycle (and can be worse during geomagnetic storms), Ion2 and Ion3 need to be considered, as they can cause range errors of a few to tens of centimetres (Wang et al., 2005).The impact of the higher order ionospheric errors on the estimated station coordinates is studied by comparing the coordinates estimated by PPP performed with the original and corrected observation files, as done in this work.
Section 1 of this paper gives an introduction of the work performed; Sect. 2 presents the literature review; Sect. 3 introduces the methodology for calculating the higher order ionospheric error terms by using the RINEX HO software and for PPP; Sect. 4 presents the results for the calculated values of Ion2 and Ion3, whereas Sect. 5 for the PPP results from both the original and corrected observation files.The paper concludes with discussion and suggestions for future work in Sect.6.

Literature review
Previous works related to the higher order ionospheric effects consider the ionospheric refractive index to derive the error due to the ionosphere and the Chapman theory for the layered structure of the ionosphere; they account for the influence of the geomagnetic field on the refractive index of the (anisotropic) ionosphere and some may neglect the differential (frequency and satellite elevation angle dependent) bending effect on the GNSS signals.Different authors consider these concepts differently to estimate the contribution of the higher order ionospheric effects to the GNSS error budget.In Wang et al. (2005) a multi-GNSS approach is taken to estimate the higher order error terms; and the authors focus on the techniques of eliminating/estimating the ionospheric errors through new linear combinations possible with the new signal frequencies of the modernized GNSS.Brunner and Gu (1991) observe that the RREs due to Ion2 and Ion3 in the dual frequency solution (i.e. using the IF observable) can reach several centimetres at low satellite elevations when the ionospheric electron density is as high as during the active period of the solar cycle.Their proposed model (using two separate Chapman functions to represent the top and bottom sides of the ionospheric electron density profiles) can eliminate the RREs to better than 1 mm by considering a series expansion of the ionospheric refractive index, an accurate ionosphere model (that provides the electron density as a function of height in the ionosphere), the International Geomagnetic Reference Field (IGRF) and also by accounting for the differential bending effect (important especially at low elevation angles) of the GPS signals (along with the tropospheric effect on the curvature of the GPS signals).Bassiri and Hajj (1993) propose an approach which can eliminate the RREs to the millimetre level by considering a series expansion of the ionospheric refractive index, a thin shell model for the ionosphere (as a sum of E, F 1 and F 2 Chapman layers), a tilted dipole model for the geomagnetic field and by neglecting the bending effect on the GPS signals (since they assume that the bending effect is insignificant for the satellite elevation angles greater than 30 • ).Ioannides and Strangeways (2002) show an analytical perturbation technique to determine the Ion2 and Ion3 terms for which they account for the ray bending effect, and the authors compare these results with those obtained from precise ray-tracing calculations for the GPS frequencies.They conclude that the refracted geometrical path increases compared with the LoS and there is a corresponding increase in the TEC with an associated phase advance.If the influence of the magnetic field for the L band signals is neglected then the total curvature error is of comparable magnitude with the increase in the geometrical path length related with the longer curved path but of opposite sign; this represents the phase advance.The authors thereafter suggest that both terms do not need to be determined since the total curvature error is of the same magnitude but opposite sign of the increase in the geometrical path.Kedar et al. (2003) focus on the impact of Ion2 on PPP by considering a co-centric tilted magnetic dipole and the GIM software (Global Ionospheric Mapping software from the Jet Propulsion Laboratory -JPL, 2010) which provides twodimensional electron density maps for the ionosphere taken as a thin layer at 400 km altitude.They use the satellite clock and orbit products not corrected for Ion2.In their comparison of the coordinate time series corrected for Ion2 with the original uncorrected coordinates, they find sub-centimeter level error contribution by the Ion2 term to the GNSS positioning error.Wang et al. (2005) present a triple-frequency method for correcting Ion2 and propose an ionosphere-free linear combination method based on three frequencies, claiming that their triple-frequency method can correct the effects to the millimetre level.Moreover, they derive a formula for Ion3 for which they apply the semi-empirical ionospheric model developed by Anderson et al. (1987) who define TEC as a linear function of the maximum electron density (N max ) in the ionosphere (Eq. 1) which gives N max as 4.405 × 10 −6 TEC; this agrees well with the linear interpolation N max ≈ 4.415 × 10 −6 TEC applied by Fritsche et al. (2005).After obtaining TEC from pseudorange (PR) measurements with L1 and L2, Wang et al. (2005) can estimate Ion3 with an accuracy of ±1 mm.Fritsche et al. (2005) investigate the impact of correcting the satellite orbits and Earth rotation parameters while estimating the station coordinates in a non-fiducial PPP approach using the Bernese GPS Software V5.0 (Bernese, 2007).Following the mathematical model of Bassiri and Hajj (1993) for Ion2 and Ion3 and using a thin shell model for the ionosphere, they apply GIMs for TEC data and a co-centric tilted magnetic dipole for the geomagnetic field.They show that both the station coordinates and the satellite orbits can change at the centimeter level when the corrections for Ion2 and Ion3 are applied.They emphasize that a consistent correction method for RREs should use the corrected GPS observations and products rather than the corrected observations without taking into account the RREs for the products.Hoque and Jakowski (2007) quantify the residual "phase" error due to Ion2 and neglect that due to Ion3 (differential bending of the GNSS signals also neglected) claiming that on a disturbed day (e.g. about 100 TECU) the RRE due to Ion3 is at the sub-millimetre level.Their model, which can provide better than 2 mm accuracy for GNSS users in Germany, does not require knowledge of the instantaneous geomagnetic field since they take the geomagnetic field component for a reference user position in central Germany.Knowledge about the electron distribution along the propagation path is also not required.These assumptions make the model applicable for real time GNSS applications in central Germany.Kim and Tinin (2007) use perturbation theory to study the residual error in the dual frequency ionosphere free observable; they investigate in particular the Ion3 term associated with ray bending effect on the GNSS signals penetrating through an inhomogeneous ionosphere.Taking into account that Ion3 term includes not only the quadratic correction due to the refractive index but also the correction for the ray bending effect, they show that the ray bending effect may dominate the Ion2 error contribution.They consider both the regular large scale and random small scale irregularities in the ionosphere such that the latter can, at times, cause residual error comparable to or greater than that of the Ion2 term, dominating the contribution to the residual error in the IF observable.Pajares et al. (2007) consider the impact of Ion2 on the satellites clocks and show that the estimates of the RREs on the receiver coordinates, satellite positions and clocks are correlated.Regarding the receiver positions, they infer that Ion2 has a sub-daily impact of less than 1 mm during March in 2001 -a year during the peak of the solar cycle.As for the satellite positions, they show that Ion2 causes a daily mean global southward displacement of several millimetres, depending on the ionization level in the ionosphere.Regarding the satellite clocks, which are most affected by the higher order ionospheric effects according to their results, RREs can cause deviations even larger than 30 picoseconds (corresponding to about 1 cm) depending on the latitude and local time of the receiver position.Datta-Brua et al. (2008) show that, unlike the Ion1 term which has the same magnitude but opposite signs for the group and carrier phase measurements, the Ion2 and Ion3 have different magnitudes and signs for these two types of measurements.For this reason, the authors claim that the higher order errors accumulate in the carrier smoothing of the IF (to the first order) code observable; they authors show that the errors in the carrier-smoothed code measurements are mostly due to Ion2 and Ion3.In other words the unaccounted higher order group errors contribute to the error in the carrier smoothing.Although can be neglected in many applications, these residual errors can be crucial in high precision applications.Pajares et al. (2008b) focus on Ion2 and different methods to obtain slant TEC (STEC) and the geomagnetic field com-ponent projected onto the receiver-satellite path.Considering the error due to Ion2 on the positioning, they emphasize that correction for Ion2 must be applied to all fiducial coordinates -application only to the unknown station (user) can lead to errors in the estimated coordinates that can be worse than applying no correction at all at the any receiver involved.Kim and Tinin (2011), in a more recent study, explore the possible ways of eliminating the higher order ionospheric error terms considering a multi-frequency GNSS approach and show how the GNSS accuracy can be improved considering the propagation of the signals through an inhomogeneous and random structure of the ionosphere.Through numerical simulation they show that the systematic, residual ionospheric error can be significantly reduced (under certain ionospheric conditions) through triple frequency GNSS.Moore and Morton (2011) focus on the Ion2 term introduced by the interaction between the GNSS signal and the magnetic field of the Earth.The anisotropy of the ionosphere causes the right hand circularly polarized (RHCP) GPS signals propagate in two (ordinary and extraordinary) modes, as a linear combination of them, depending on the angle between the GPS signal wave normal and the Earth's magnetic field.These two modes correspond to two different magneto-ionic polarizations each with a particular refractive index that needs to be considered in the Ion2 term.The authors show that near the geomagnetic equator signals arriving from the north propagate with the ordinary polarization (associated with left hand circularly polarized wave) yielding a "positive" Ion2 term for the carrier phase; and those arriving from the south propagate in the extraordinary mode polarization (associated with right-hand circularly polarized wave) making the Ion2 term "negative" for the carrier phase.A positive Ion2 term corresponds to the presence of error still to be accounted for in the (first order) ionospherically corrected measurements.The authors also point out a misconception in the work of Bassiri and Hajj (1993) who assume that the left hand circularly polarized (LHCP) component of GPS signals propagates in the ordinary mode and do not realize the fact that the RHCP signal component may indeed travel in either of the propagation modes.Moore and Morton (2011) show that the magneto-ionic polarization of the predominantly RHCP GPS signal depends on the direction of the GPS signal wave vector with respect to the magnetic field line.Considering three different geographic locations to show the influence of this fact on the Ion2 term, Moore and Morton (2011) show that Ion2 is asymmetric with respect to the geomagnetic equator such that depending on the magnitude of the angle between the wave vector and the magnetic field line, a RHCP wave has different propagation modes thus considering only one propagation mode is expected to lead to mismodelling inaccuracies in estimating the error due to the Ion2 term.
From the literature review it can be understood that, while estimating the magnitudes of the errors due to the Ion2 and Ion3, higher accuracy can be achieved by (1) using a more precise geomagnetic field like the IGRM instead of the dipole model; (2) obtaining accurate estimates for the electron content along the signal link (STEC values) which can be either retrieved from Global Ionospheric Maps (GIMs) or estimated from PR measurements; (3) using the satellite and orbit products estimated by applying corrections for the higher order ionospheric error terms (this is particularly important for a systematic and accurate analysis of the impact of the higher order terms on PPP).Regarding point 2, vertical TEC (VTEC) data from GIMs can be converted into STEC by making use of a single layer mapping function in RINEX HO.Alternatively STEC can also be obtained from the PR measurements on the L1 and L2 frequencies (Eq.2); this, however, requires inputting the receiver and satellite Differential Code Biases (DCB rec and DCB sat , respectively) which were obtained from CODE for use in RINEX HO in this work.The uncertainty in any of the terms on the right hand side of (Eq.2) propagates into the calculation of STEC on a particular receiver-satellite link, according to the error propagation law (Eq.3).Although STEC can be estimated with a comparable accuracy using either GIMs or PRs (Marques et al., 2007), the non-availability of the DCB values at some instances hinders estimation of STEC from PRs. Therefore STEC is obtained from GIMs in RINEX HO in this work.Regarding the third point, since such corrected products were not available during the progress of this work, the satellite and orbit products estimated without corrections for the higher order error terms were used.

Methodology
The observation stations in this work are selected from the International GNSS Service (IGS, 2010) network, aiming for a reasonably good latitudinal (mid and high latitudes, including the auroral region) and longitudinal coverage in Europe (Fig. 1); the stations coordinates are provided in Table 1.Four sets of days (Table 2) are selected for analysis.In order to investigate the impact of the solar activity devoid of disturbances in the geomagnetic field, day-of-year (DOY) 312-316 in 2001 and 321-326 in 2006 were selected; for these periods, the planetary geomagnetic index, K p , is <4, which is a good threshold to exclude the influence of geomagnetic storms (NOAA, 2010).In order to include the influence of a more disturbed geomagnetic field, DOY 294-296 in 2001 and 301-307 in 2003 were selected, when K p was ≥4.Other geomagnetic indices like the AE (Auroral Electrojet index) or Dst (Disturbance Storm Time index) could also be considered (World Data Center for Geomagnetism -WDC, 2011) while selecting the days for analysis.However the K p index was deemed adequate, given the focus of this work on the mid-to-high latitudes.
The range errors in the GPS observables due to Ion2 and Ion3 (on GPS L1 and L2 frequencies) were estimated and corrected using the software tool RINEX HO (Marques et al., 2011), developed at the Sao Paulo State University in Presidente Prudente, Brazil.The program requires as input the observation and navigation files (in the receiver independent exchange, RINEX, format), and GIM files or DCB information (according to the user's choice of the method for TEC calculation).An input text file is used, with the relevant file names and execution specifications (in this case the choice of the method to calculate TEC).The program applies the corrections to the GPS code and phase observables, corrects the input observation file accordingly and returns the output files (corrected observation file, corrections for Ion2 and Ion3 on L1 and L2 frequencies for code observations and STEC for each receiver-satellite link).
The corrected observation file allows to perform positioning in order to compare the station coordinates estimated using the corrected and original observation files.This allows assessment of the impact of the higher order ionospheric effects on PPP, which is accomplished using the Bernese V5.0 (Bernese, 2007) software in this work.2 2 1 1 2 1 0 1 18 November 322 0 0 0 0 1 0 0 1 19 November 323 1 1 1 1 1 0 0 0 20 November 324 0 0 0 0 1 1 0 1 21 November 325 0 0 0 0 0 0 0 1 22 November 326 0 0 0 0 1 1 3 3 The code and the carrier phase GNSS observation equations are given in Eqs. ( 4) and (5), respectively.The ionospheric delay term (I f 1 ) appears with a "+" sign for the code delay (Eq.4) and with a "−" sign for the advance of the carrier phase (Eq.5).Focusing on the PR measurements (Eq.4), the ionospheric delay effect (I f 1 ) can be represented more explicitly as a series in inverse powers of frequency within the geometric optics (GO) approach i.e. considering the refraction of the GNSS signals penetrating through (large scale) electron density irregularities in the ionosphere.A series representation of the delay effect on the PR measurements, δρ g in Eq. ( 6), can be derived when the Appleton-Hartree equation (Budden, 1966) is considered for the ionospheric refractive index that is different than unity.Based on the GO approach, the Appleton-Hartree formula leads to the three terms on the right hand side of Eq. ( 7) which are the first (Ion1), second (Ion2) and third (Ion3) order ionospheric effects that delay the code measurement, respectively.The same order effects for the carrier phase measurements are given in Eq. ( 8).Hereafter the discussion considers the ionospheric effects on the PR measurements (Eq.7); the arguments are however applicable for the carrier phase range measurements taking into account the correct sign notation and coefficients for these three terms.N ds term in Eq. ( 6) is the integral of the electron density (N) along LoS between the receiver and satellite inside a columnar cylinder of unit cross sectional area such that the integration gives the (slant) total electron content, along LoS, STEC (1 TEC unit, 1 TECU, is 10 16 e − /m 2 .At times of low background solar activity, as during the quiet periods of the solar cycle, STEC is usually around 20-30 TECU at mid latitudes, corresponding to about 3-8 m range delay on GPS L1 frequency giving negligible RREs due to Ion2 and Ion3 terms under these conditions (1 TECU has about 0.16 m delay effect on GPS L1, Kintner Jr., 2006).STEC depends on the geometry of the receiver and satellite link, time of day (with a diurnal variation that attains a peak around local noon), time of year (seasonal dependency) and the solar cycle (greater STEC during peak of the solar cycle).The diurnal variation of STEC on the receiver-satellite links for GPS L1 can be seen for the observation stations in Figs. 2, 5, 8 and 11, for the different levels of solar and geomagnetic conditions.B 0 cosθ term in the Ion2 (Eq. 6) can be taken out of the integral assuming that it is LoS-independent, leaving N ds, i.e.STEC.Ion3 term contains the integral of the square of the electron density, N 2 ds (Eq.6), which can be estimated using the shape parameter η (Hartmann and Leitinger, 1984).This helps to approximate the ionospheric electron density profile in terms of the maximum electron density, N max , and the shape parameter, η, giving ηN max STEC for this integral.The shape parameter η can be taken as 0.66 which is valid for different satellite elevations and maximum electron densities (Hartmann and Leitinger, 1984).Based on these arguments, Eq. ( 6) can be written as Eq. ( 7).For the carrier phase measurements, Eq. ( 8) represents the ionospheric range errors (in meters) up to the third order, neglecting the bending effect that is associated with the third order error term   Eq. ( 7), the magnitudes of Ion2 and Ion3 in the code-based measurements can be expressed as in Eq. ( 9) and Eq. ( 10), respectively.From Eqs. ( 7) and ( 8), it can be seen that the higher order ionospheric error terms for code measurements can be estimated from the carrier phase measurements, and vice versa, applying appropriate multiplicative terms (e.g.magnitude of Ion2 for code observations is twice as large as that for carrier phase observations) and sign notation for each order term.Due to the multiplicative terms it can be understood that the first order linear combination of PR observations does not eliminate the higher order error terms.It can also be seen (Eqs.7 and 8) that TEC along LoS is important for all the error terms (it should be reminded that the bending effect in Ion3 is neglected here).Moreover, due to the LoS dependency of TEC the magnitudes of Ion2 and Ion3 change according to the receiver-satellite geometry.These residual error terms become more important in the differential positioning mode (especially for long baselines when signal links pierce through comparably different parts of the ionosphere) and at low latitudes, in particular during high solar activity (when ionization in the ionosphere is expected to be greater).
Accurate values of STEC for Ion2 and Ion3 can be obtained from (a) Global Ionospheric Maps, GIMs, which contain VTEC data accurate to about 2-8 TECU (Feltens and Schaer, 1998) or (b) PR measurements according to Eq. ( 2).These two methods show comparable accuracy (2-8 TECU) in the estimated STEC values (Marques et al., 2007).For the former method, VTEC from GIMs is converted into STEC using a mapping function, considering the IPP, given by the intersection of the receiver satellite path with the ionosphere assumed as a single thin shell at an altitude of 450 km (same value as taken in RINEX HO) according to the receiver-  satellite LoS geometry.As stated by Pajares et al. (2008b), GIMs can provide less accurate STEC values at the low latitudes for low elevation satellites; yet, since in this work midlatitude stations and a cutoff angle of 10 • are considered, this is not of concern here.For the latter, there is the need for the receiver and satellite interfrequency biases (also referred to as Differential Code Biases, DCB rec and DCB sat , respectively).These frequency-dependent biases are relatively constant in time and must be input in RINEX HO.Within this work they were provided from the Center for Orbit Determination in Europe (CODE, 2010).Non-availability of these biases may halt the process of STEC estimation from PRs (Eq.2).Thus for continuity of calculations STEC values were obtained from GIMs (a user option in the program).
As can be seen in Eq. ( 9), Ion2 depends on the projection of the geomagnetic field (B) onto the receiver-satellite link, (B 0 cosθ), which can be more accurately calculated if a precise geomagnetic field like the International Geomagnetic Reference Model (IGRM, 2010) is used instead of a dipole model.The GEOPACK library (Geopack subroutines, 2011) contains FORTRAN subroutines for computing the geomagnetic field in the Earth's magnetosphere, transforming between various coordinate systems and tracing along field lines (Tsyganenko, 2001).IGRM is used in RINEX HO for a physically more realistic and accurate modelling of the geomagnetic field.
In Eq. ( 10) it can be seen that Ion3 depends on N max for which a linear interpolation can be used to approximate N max in terms of TEC (Fritsche et al., 2005).A modified version of this interpolation is also available (Piraux et al., 2010) where N max is redefined in terms of TEC (Eq. 11).
The RREs (due to Ion2 and Ion3) are calculated for GPS L1 and L2 carrier frequencies in RINEX HO; however the  Based on the equations for the higher order range errors (Eqs.7 and 8), a straightforward calculation (taking 150 TECU along line of sight, B 0 cosθ about 2.7 × 10 −5 Tesla at IPP taking B 0 as 3.12 × 10 −5 and θ as 30 • and N max as 4.416 × 10 −6 TEC) gives about 24 m, 2.3 cm and sub-millimeter (negligible which may be due to the fact that bending effect is excluded in this calculation) level range errors for Ion1, Ion2 and Ion3, respectively, for GPS L1.For such conditions, for instance, Ion2 is about 0.09 % of Ion1 for the code based range measurements (equivalently 0.05 % for carrier phase range measurements).In this case, using the IF observable to remove the Ion1 term would be accurate up to about 99.9 %; this agrees well with the estimation of Klobuchar (1987) that the IF observable is accurate up to about 0.1 %.However the crude values assigned to the parameters involved in Eqs. ( 7) and ( 8) should be kept in mind in this estimation.
The final stage of the work presented here is to analyze the estimated station coordinates when processing the data in PPP, using the Bernese software (Bernese, 2007).This part of the work focuses on the impact of using the corrected (for Ion2 and Ion3) observation files in PPP in order to infer the significance of the contribution from the higher order ionospheric effects in GNSS positioning in the European region during different background physical (solar, geomagnetic, ionospheric) conditions.
In the coordinates estimation process, Ion2 and Ion3 corrections are applied only to the receiver observations; the orbit and clock products (used in the Bernese software) in-  volved in PPP are computed from a global network that does not apply corrections for Ion2 and Ion3.It must be noted that a systematic and accurate investigation of the impact of the higher order terms on PPP requires the use of satellite and orbit products estimated while accounting for these higher order error terms.As shown by Pajares et al. (2008b), a more consistent and correct approach to consider (e.g.Ion2 in PPP) can be to perform dual frequency, carrier phase differential positioning where both the orbits (and other satellite products) and user coordinates are estimated considering the Ion2 correction.According to Pajares et al. (2007) the effect of Ion2 on the satellite clocks can be larger than 30 picoseconds (1 cm in range equivalent units) and several millimetres on the satellite positions.different analysis centres which may consider either Ion2 or Ion3 or both.For instance, as of 2009 JPL has been considering only the Ion2 in their ionospheric model, whereas CODE considers both Ion2 and Ion3 for their products.As stated by Pajares et al. (2008a), using the standard products, which are not corrected for the higher order ionospheric effects, with the corrected GPS observations, blurs the net impact of corrections.However, since a set of standard satellite orbit and clock products are yet unavailable (Piraux et al., 2010), the JPL satellite and orbit products estimated without accounting for Ion2 and Ion3 were used in PPP in this work, where only the receiver positions are estimated based on data corrected for the higher order ionospheric effects.

Results for the higher order error terms and the STEC values
The results (RREs for Ion2 and Ion3, PPP station coordinate differences) are shown for Ion2 and Ion3 on the code observations for the GPS L1 frequency since this applies to a wider user community using the civil code on GPS L1.Due to the inverse frequency dependency of both Ion2 and Ion3, the calculated Ion2 values for GPS L2 are about 2.11 times and the Ion3 values for GPS L2 are about 2.71 times those obtained for the GPS L1.

DOY 312-316 in 2001
These results are presented in Figs.electrons and ions tend to recombine reducing the amount of ionization, at the high latitudes night time ionization can continue due to the movement of the ionization from the daytime to the night-time part of the Earth as well as due to energetic particles arriving with the solar winds to the vicinity of the Earth and penetrating into lower altitudes of the ionosphere along the almost vertical geomagnetic field lines at these high latitudes.In the latter case, the particles moving vertically downward collide with the ionospheric particles causing collisional ionization.Such effect can be observed especially at the high latitudes where the geomagnetic field lines can route the particles and during the active period of the solar cycle when the solar radiation is stronger (Buonsanto, 1999).
As seen in Eq. ( 9) Ion2 has a LoS dependency; thus for stations at different latitudes the values for Ion2 (Fig. 3) change from being more confined to about −1 cm (for TRO1) to scattering between ±2 cm (for MATE).The negative values for Ion2 (in all plots) are due to the B 0 cosθ term, which can attain positive or negative values depending on the satellitereceiver geometry.A mid-latitude station (e.g.MATE) can track the satellites with a wider range of elevation angles whereas a high latitude station as TRO1 tracks with a more confined range of elevation angles; thus the LoS dependency and thereafter the Ion2 errors are different for the receivers at different latitudes.
In Fig. 4, it can be seen that the estimated Ion3 is significant during the noon-hours for all stations, ranging from about 2 to 3 mm from north (TRO1) to south (MATE).It can be observed that for the high latitudes, there can be significant correction for Ion3 at the night-time hours (see station TRO1 in Fig. 4).

DOY 294-296 in 2001
These results are presented in Figs. 5 to 7. Comparing the results in Fig. 5 with those in Fig. 2, it can be noticed that the night time enhancement in TEC values can be as large as the noon-time values (e.g.station TRO1 in Fig. 5) when there are geomagnetic storms in addition to high background solar activity.This causes almost two peaks for the diurnal TEC values, especially for the high latitude stations.Movement of the ionization from the day side to the night side of the Earth at the high latitudes can be enhanced by the geomagnetic storms (Ho et al., 1997).As seen in Fig. 5 for ONSA, enhancement in the auroral TEC at the night time can be due to the expansion of the auroral oval by the influence of the geomagnetic storms during this peak year of the solar cycle as also evident in the high K p values.Such night-time enhancement in TEC is not apparent in Fig. 2 for ONSA.The mid-latitude stations are observed not to have enhanced levels of night time TEC, which means that the equatorward expansion of the auroral oval was not significant enough as to influence the ionosphere above these stations during these geomagnetic conditions.
It can be seen in Fig. 6 that the error due to Ion2 is overall within 1-2 cm, which agrees to those observed in Fig. 3.However, due to the geomagnetic storms -which may cause enhancements in the ionization levels (i.e. higher TEC values), Ion2 values for ONSA can be noticeably large during the night-time hours as well -such enhancement is not apparent in Fig. 3 for ONSA.Also, the night-time values for Ion2 at TRO1 can be as large as the noon-time values (TRO1 and ONSA in Fig. 6).Since the B 0 cosθ term in Ion2 calculated by IGRM does not in general consider the actual geomagnetic disturbances, the enhanced values of Ion2 can be more correctly related with the enhancement in TEC.However it should be noted that TEC can increase due to the geomagnetic storms during such conditions.The similarity between the Ion2 error profile (Fig. 6) and the corresponding TEC values along the signal links (Fig. 5) suggests a strong dependency of Ion2 on STEC.
Compared with Fig. 4, it can be seen in Fig. 7 that when geomagnetic disturbances are also considered the error due to Ion3 becomes significant during the midday hours for all stations, however this time the magnitude of the error ranges from about 1 to 5 mm from north (TRO1) to south (MATE).The noticeable night time peaks in Fig. 7 for TRO1 and ONSA suggest that for the high latitude stations, Ion3 may need to be considered during night-time hours as well, Overall it can be seen in both Figs. 4 and 7 that Ion3 is important during daytime, amounting to a few millimetres in the range errors at these stations during the peak of the solar cycle.

DOY 301-307 in 2003
These results are presented in Figs. 8 to 10.Although within a post-peak year of solar cycle 23, the period DOY 301-307 in 2003 coincides with the so-called Halloween Storm, during which the solar radio flux index, F 10.7 (measure of radio emission from the sun at 10.7 cm wavelength correlating well with the sunspot number and used as an indicator of solar activity (Ionospheric Prediction Service -IPS, 2011)), went up to as high as 270/275 solar flux units (Space Weather Prediction Center -SWPC, 2011).It should be mentioned, however, that the change in the geomagnetic field inductance due to geomagnetic disturbances and the estimation of its influence on the Ion2 term is not straightforward.
During this period, increase in the Ion2 error term is expected to be due to higher levels of solar activity as well as to the disturbances in the geomagnetic field.Similarly, increase in Ion3 can be associated with higher levels of solar activity during this period.Results show high levels of ionization  During this period Ion2 ranges from about 1 to 2 cm (Fig. 9) and Ion3 from sub-millimetre to about 2 mm (Fig. 10) from TRO1 to MATE in both cases.Thus, even slightly outside the highest phase of the solar cycle, the solar activity can be strong enough to enhance the ionization in the ionosphere.During the high phase of the geomagnetic storm (DOY 301-302, 2003), Ion2 and Ion3 have significant enhancement; whereas during the absence of such storms (e.g.DOY 312-316, 2001) the higher order error terms are more predictable.

DOY 321-326 in 2006
These results are presented in Figs.11 to 13.In 2006, it can be seen that during the quiet period of the solar cycle, ionization in the ionosphere is low, about 20-40 TECU (Fig. 11) and the corresponding higher order range errors are also less significant than those during the active period of the solar cycle: Ion2 is observed to be one order of magnitude smaller (at millimetre level as seen in Fig. 12 as opposed to centimetre level as observed during other analysis periods) and Ion3 is at sub-millimetre level (Fig. 13) during this quiet period of the solar cycle in absence of geomagnetic disturbances.
Between the high and low periods of the solar cycle (November 2001 and November 2006, respectively, without the influence of geomagnetic disturbances in both cases) there is a significant difference in the calculated STEC values (e.g. as high as about 160 TECU during the peak of the solar cycle in November 2001 at the mid-latitudes, and about 40 TECU during the quiet period in November 2006 at the same latitudes).Significant differences are also observed in the Ion2 and Ion3 values -they change by an order of magnitude between the two periods, therefore indicating that the significantly different TEC values along the signal paths cause the different values for the higher order range errors during these two periods.This highlights the importance of considering the higher order range errors during the upcoming solar maximum, predicted for around 2013.
During this quiet period of the solar cycle (i.e.characterised by low levels of ionization in the ionosphere), when geomagnetic field disturbances are negligible the high order ionospheric error terms should not degrade the measurement accuracy significantly.

PPP results
PPP is a high accuracy positioning method which can be performed with a dual frequency receiver (so that the IF observable can be used) and that exploits the use of highly accurate externally provided (e.g. by the IGS) satellite orbit and clock corrections (Bernese, 2007).Due to its high (potentially centimetre level) positioning accuracy, PPP was performed in this work with the Bernese software to investigate the impact of correcting the GPS range observations for the errors due Ion2 and Ion3.With the Bernese software, a free network solution can be carried out, i.e. a solution where the satellites orbits define the coordinates system to which the estimated positions refer to.
As detailed before, RINEX HO applies corrections for Ion2 and Ion3 to the GPS observation files in the RINEX format and outputs a corresponding corrected observation file in the same format.In this work, PPP was performed respectively with both observation files, i.e. a file with the "uncorrected" and another with the "corrected" (for the higher order terms) observations.Figures 14 to 17 show how much the stations coordinates differ in both cases (in latitude, longitude and ellipsoidal height) for all four sets of days and stations being analysed.The differences in all three coordinate components are computed by subtracting the PPP results obtained with the original observation (uncorrected) files from those obtained with the corrected files (see Table 3 for the numerical values of the differences).The results are discussed below: -Considering the high solar activity period with negligible disturbances in the geomagnetic field (DOY 312-316 in 2001, Fig. 14) when the corrections for Ion2 and Ion3 account majorly for the impact of the solar activity, it can be observed (Fig. 14, top plot) that the high latitude stations get northward corrections (about 2-3 mm) and mid-latitude stations southward (about 1 cm).During this period all stations are observed to have westward corrections (about 1-2 cm) in general (Fig. 14, middle plot) and the vertical component of the station coordinates were greater (by about 2-3 cm) in general for the mid-latitudes and smaller for the high latitudes (Fig. 14, bottom plot) when the corrected observation files were used in PPP.Pajares et al. (2007) who focus only on the Ion2 term and its impact on the geodetic estimates show that applying Ion2 correction to subdaily differential positioning (using IGS data network) changes the receiver positions at submillimeter level which is northward for the high latitudes and southward for the low latitudes.
-Considering the high solar activity period with disturbed geomagnetic conditions when K p was as large as 7 (DOY 294-296 in 2001, Fig. 15), it is difficult to decide on a general common trend for the latitudinal corrections (of few millimetres, see Fig. 15, top plot) in general (Fig. 15, top plot).During this period, the previously (Fig. 14, middle plot) observed westward correction seems to be suppressed.The estimated vertical corrections (Fig. 15, bottom plot) are overall upward however the mid-latitudes are corrected downward (at sub-cm level) on average.It should be pointed out that the short observation period considered here may hinder a more conclusive analysis.
-During the post-peak period of the solar cycle with disturbances in the geomagnetic field, during the so-called Halloween Storm, (DOY 301-307 in 2003, Fig. 16), overall a southward correction (Fig. 16, top plot) can be deduced with magnitudes mostly at millimetre level but at times a few centimetres.A distinguishable feature during this post-peak period can be observed in the longitudinal corrections (Fig. 16, middle plot): the geomagnetic activity seems to suppress changes in the longitude component of the station coordinates.Regarding the estimated heights of the stations during this period, 1-2 cm level corrections can be observed (Fig. 16, bottom plot) such that the high latitudes are corrected downward and mid-latitudes upward.
-Considering the period of low background solar activity with quiet geomagnetic conditions (DOY 321-326 in 2006, Fig. 17), it can be seen in the horizontal station components that PPP results do not show significant differences when the observation files are corrected for Ion2 and Ion3 (Fig. 17, top and middle plots).It is difficult to observe a general trend in the direction for the vertical corrections (Fig. 17, bottom plot).
Considering that PPP can potentially provide centimetre level accuracy for the estimated station coordinates and that the corrections for Ion2 and Ion3 are about centimetre and millimetre levels, respectively, during adverse ionospheric and geomagnetic conditions, it can be expected that impact of Ion3 in PPP may be unnoticeable due to the noise level of the positioning.It can be expected that the differences in the estimated station coordinates (using the corrected and original observation files) would be mostly due to corrections for Ion2 and then for Ion3 and the noise level of the positioning solution.
6 Discussion and suggestions for future work Based on the analysis involved in this work the following points can be drawn: -Enhancement in TEC values can occur due to greater solar activity, as during the peak of the solar cycle.Similarly the geomagnetic field disturbances can also drive mechanisms that enhance ionization in the ionosphere.As TEC is an important parameter while calculating the range errors due to the ionosphere, more contribution from the ionospheric error terms is expected when the TEC values are higher.During the post-peak period of the solar cycle, if geomagnetic field disturbances are present (e.g.DOY 301-307 in 2003), the higher order ionospheric error terms were observed to contribute to the overall measurement accuracy at magnitudes comparable to those occurring during the peak of the solar cycle without the presence of such disturbances (e.g.DOY 312-316 in 2001).During the quiet period of the solar cycle (low ionization levels in the ionosphere), when the geomagnetic field disturbances were negligible (e.g.DOY 321-326 in 2006), the higher order ionospheric error terms were observed to be very small.Even during the post-peak years of the solar cycle, high TEC values may be observed which can be explained by the geomagnetic activity that can route the incoming solar particles to lower altitudes in the ionosphere, especially at the high latitudes where the geomagnetic field lines are mostly directed towards the surface of the Earth.Thus, enhancement in TEC should not be expected only during the peak years of the solar cycles; it can be observed in general correlated with increased levels of geomagnetic activity.
-In terms of the effects of the Ion2 and Ion3 in longitude, a general westward correction was observed in this work (during active period of the solar cycle) where the mid-latitude stations were observed to be affected more than the high latitude ones  and northward for the high latitudes.For the height component, the mid latitude stations are observed to be corrected upward in general and the high latitude ones downward.
-Diurnal variation in RREs was such that a minimum was observed before sunrise and after sunset, with a maximum around noon for the stations analysed in Europe.
The strong diurnal variation in TEC is expected to be a reason for this.
-Had the corrected satellite orbit and clock products been used in PPP, a more systematic and realistic analysis of the differences in the positioning results could have been carried out.In the approach followed in this work, the net effect of correcting the observation files is expected to be obscured since the corrections were applied only to the receiver observations.
-Comparing the Ion2 and Ion3 values presented in this work, the diurnal variation in Ion3 is stronger than that in Ion2 (i.e. the relative difference between the minimum and maximum for Ion3 and greater than that for Ion2).This can be explained by the fact that in addition to the dependence on STEC, Ion2 depends on the projection of the geomagnetic field onto the signal path and Ion3 on the maximum electron density in the ionosphere; in the former no diurnal variation is anticipated (IGRM consideration of TEC for the magnetic field is not there thus a clear relation between this magnetic term and TEC cannot be established), whereas in the latter diurnal variation can be anticipated since the maximum electron density normally reaches a maximum at about the local afternoon (Ratcliffe, 1972).-It was observed that both the horizontal and vertical components of the estimated station coordinates in PPP can differ at cm-mm level when the corrected (for Ion2 and Ion3) observation files are used.Differences in the estimated horizontal components can be influenced by not only the presence but also the strength of the geomagnetic field disturbances.It was also observed that when PPP is performed respectively with both the corrected and original observation files, the differences in the estimated station coordinates during the non-peak period of the solar cycle can be comparable to those during the peak period if the former has contribution from significant levels of geomagnetic field disturbances.
-This work describes a method to analyze the higher order ionospheric effects on the GNSS observations and on positioning; future work will be carried out on the basis of this methodology, further accounting for the bending effect in the Ion3 term.Longer term periods will also be analyzed -monthly or yearly analyses for the higher order ionospheric effects can be performed with open-sky observations during the upcoming solar maximum, expected around 2013.Future work can also consider the new GNSS signals like GPS L2C, L5 and Galileo L1 and E5 which may be available by then.This can allow an assessment of the advantage of these new signals for the GNSS user community to correct range errors related to the ionosphere, especially its higher order terms that have so far been mostly ignored.
Fig. 16."Corrected -Uncorrected" PPP results for the geodetic coordinates for the stations on DOY 301-307 in 2003.
Fig. 16."Corrected -Uncorrected" PPP results for the geodetic coordinates for the stations on DOY 301-307 in 2003.

Table 1 .
Coordinates of the IGS stations considered for analyses in this work.

Table 2 .
Days used in the analyses (given as day-of-year, DOY) and the corresponding 3-hourly K p values for each day such that the first K p value corresponds to midnight 00:00 LT.
These products can be obtained from Ann. Geophys., 29, 1383-1399, 2011 www.ann-geophys.net/29/1383/2011/ 2 to 4. The diurnal variation in STEC values can be seen in Fig. 2 for all stations.The midday peak values are greater at the mid-latitudes (e.g. at MATE) than at the high latitudes (e.g.TRO1).For the high latitude station TRO1, the night-time enhancement in TEC values can be as large as half the noon-time values.Although ionization due to solar radiation is absent at night and free

Table 3 .
Differences in the calculated station coordinates (delta height/latitude/ longitude, in meters) when PPP is performed with the corrected and uncorrected observation files.Positive differences in height, latitude and longitude are upward, northward and eastward, respectively.