Simulation study of the inﬂuence of the ionospheric layer height in the thin layer ionospheric model

This work aims to contribute to the understanding of the inﬂuence of the ionospheric layer height (ILH) on the thin layer ionospheric model (TLIM) used to retrieve ionospheric information from the GNSS observations. Particular attention is paid to the errors caused on the estimation of the vertical total electron content ( vTEC ) and the GNSS satellites and receivers inter-frequency biases (IFB), by the use of an inappropriate ILH. The work relies upon numerical simulations performed with an empirical model of the Earth’s ionosphere: the model is used to create realistic but controlled ionospheric scenarios and the errors are evaluated after recovering those scenarios with the TLIM. The error assessment is performed in the Central and the northern part of the South American continents, a region where large errors are expected due to the combined actions of the Apple-ton Anomaly of the ionosphere and the South-Atlantic anomaly of the geomagnetic ﬁeld. According to this study, there does not exist a unique ILH that cancels the vTEC error for the whole region under consideration. The ILH that cancels the regional mean vTEC error varies with the solar activity and season. The latitude-dependent conversion error propagates to the parameters of the model used to represent the latitudinal variation on the vTEC on the ionospheric layer, and to the IFB, when these values are simultaneously estimated from the observed sTEC . Besides, the ILH that cancels the regional mean vTEC error is different from the one that


Introduction
The total electron content (T EC) extracted from the GNSS measurements is a valuable source of information for a great variety of ionospheric studies (e.g., Hernández-Pajarez et al. 2009).T EC is defined as the curvilinear integral of the electron density distribution, N , along a given trajectory, (e.g., Davies and Hartmann 1997): where N is measured in electron per cubic metre, in metre, and T EC in total electron content unity (TECu).The notation vT EC is preferred when is a vertical line, while sT EC is preferred when is a slant line.
In order to provide absolute T EC values, GNSS observations must be calibrated to remove the inter-frequency biases (IFB) originated in the satellite and the receiver hardware and firmware (e.g., Sardon et al. 1994).In most cases this calibration is done using the so-called thin layer ionospheric model (TLIM), in which the whole ionosphere is represented by one spherical layer of infinitesimal thickness with equivalent vT EC, situated somewhat above the F2 ionospheric peak, between 350 and 450 km above the Earth's surface (e.g., Mannucci et al. 1998).Within this model, the sT EC along a given line-of-sight (LOS) from a GNSS satellite, S, to a GNSS receiver, R (Fig. 2), and the vT EC L at the point where the LOS crosses the layer (the so-called ionospheric penetration point, IPP), are related by: where z L is the LOS zenith angle at the IPP and sec z L is the mapping function (MF) used to convert the slant into vertical T EC.Furthermore, the vT EC L is parameterized as a function of time, t, and the geographic latitude, ϕ L , and longitude, λ L , of the IPP: vT EC L = f L (t, ϕ L , λ L ; x 0 , . . ., x n ), where x i , i = 0, . . ., n, are the function parameters that are estimated from the GNSS observations based on the following equation of observation: where L I is the dual-frequency GNSS ionospheric observable, the so-called geometry-free linear combination that is formed by subtracting simultaneous observations at the two different frequencies L1 and L2, and in this way removing all frequency-independent effects, v I is the associated observational error, and β = β R + β S is the sum of the satellite and the receiver IFB that should be estimated from the observations together with the function parameters.
The ionospheric layer height above the Earth's surface, h L , is the key parameter for the TLIM (e.g., Klobuchar 1996), it defines the MF: where z E is the LOS zenith angle at the Earth's surface, r L = r E + h L , and r E is the mean Earth's radius; besides, h L defines the geographic coordinates of the IPP ϕ L and λ L .
In spite of its simplicity, the TLIM has been and still is widely used in a variety of problems such as: correcting single-frequency GPS observations using the ionospheric information broadcasted by the GPS satellites in connection to the Ionospheric Correction Algorithm (often called the Klobuchar's model) (Klobuchar 1987); computing the ionospheric correction in the USA WAAS (Wide Area Augmentation System; e.g., Enge et al. (1996); computing the International GNSS Service (IGS) Global Ionospheric Maps (GIM; e.g., Schaer 1999), etc.
Since the choice of the ILH has a direct impact on the computed vTEC L value, the problem of determining the ILH for the TLIM has been studied and a range of different ILH can be found in the literature (e.g., Mannucci et al. 1993;Ciraolo and Spalla 1997;Birch et al. 2002).This paper intends to complement previous efforts (Komjathy and Langley 1996; Radicella et al. 2004;Nava et al. 2007;Smith et al. 2008;Brunini and Azpilicueta 2010) to characterize the influence of the ILH on the TLIM and to assess their effects on the vTEC L and IFB estimation.The focus of this work is to revisit the problem and to provide a detailed assessment of the magnitude of the effect due to the use of the TLIM.For studying how this effect varies according to solar activity and season, the assessment was done on a monthly base and covering from low to high solar activity trying to represent all possible real scenarios for the ionosphere in one of the most affected regions of the planet.The simulation study is based on the use of an empirical model of the Earth's ionosphere, which allows creating realistic but controlled ionospheric scenarios for the evaluation of the errors that are produced when the TLIM is used to reproduce those scenarios.The error assessment is performed for the Central American and the northern part of the South American continents, where largest errors are expected because the combined actions of the Appleton Anomaly of the ionosphere and the South-Atlantic anomaly of the Earth's magnetic field.

Ionospheric scenario for the error assessment
The NeQuick ionospheric model is used to create the ionospheric scenarios.The first version of NeQuick has been used to create ionospheric scenarios for the EGNOS project and is being implemented for single-frequency applications of Galileo (Radicella and Leitinger 2001).It has also been adopted by the International Telecommunication Union, Radiocommunication Sector (ITU-R) as a suitable method for total electron content (TEC) modeling (ITU 2003).The NeQuick model has been deeply studied to validate its formulation and to obtain improvements, (Coïsson et al. 2006;Leitinger et al. 2005).These studies and recent improvements originated the second version of NeQuick (Nava et al. 2008), which is used for this research.NeQuick computes the ionospheric electron density as a function of solar activity, month, UT, height and geographic coordinates.It is a quick-run model that allows to calculate vertical and slant TEC for any specified path.From 90 km up to the F2 peak, this model uses a modified DGR profile formulation (Di Giovanni and Radicella 1990) which includes five semi-Epstein layers with modelled thickness parameters.It is also based on anchor points defined by foE, foF1, foF2 and M(3000)F2 values.The model topside is represented by a semi-Epstein layer with a height-dependent thickness parameter (Hochegger et al. 2000) empirically determined (Coïsson et al. 2006).The NeQuick (Fortran 77) source code is available at: http://www.itu.int/oth/R0A04000018-/en.
In this work, the 'ground-truth' is provided by the Nequick model.No inaccuracies in the model are considered since the study is centered in the assessment of the errors due to the sTEC to vTEC conversion by means of a commonly used MF.The error assessment is performed in the Central American and the northern part of the South American continents.Due to its proximity to the Earth's Equator, these regions receive a great portion of the solar radiation, which causes an intense ionization process.The action of thermospheric winds creates a West to East electric field, E, which in combination with the almost horizontal North to South geomagnetic field, B, produces an uplift of the free electrons over the geomagnetic equator as a consequence of the E × B vertical drift.The free electrons move up and then away from the geomagnetic equator, to the South and the North, and then go down following the geomagnetic field lines.This so-called "fountain effect" creates two regions of height ionization in both sides of the geomagnetic equator and, consequently, two crests in the vT EC distribution, which are known as the Appleton (or Equatorial) Anomaly (Appleton 1946).Besides, the South Atlantic Anomaly creates a severe distortion of the geomagnetic field, which makes the morphology of the Appleton Anomaly even more complex (Abdu et al. 2005).
Since the largest horizontal gradients are expected to happen in the North-South direction at the time of maximum evolvement of the Appleton Anomaly crests, the ionospheric scenarios created for this research are latitudinal cuts of the ionosphere through the −60 • meridian, at 14 LT (Fig. 1).

Methodology for the error assessment
Let the plane of Fig. 2 be coincident with a geographic meridian; let G be the geocentre and G Q the projection of the Earth's Equator plane on that meridian plane; and let's assume a hypothetical GNSS receiver, R, located just on that meridian plane at geographic latitude ϕ R , and a hypothetical GNSS satellite, S, located just on the same meridian plane.The position of any point in that plane can be defined by the polar coordinates, r = r E + h and ϕ, and the coordinates of any point, P, along the GNSS satellite-to-receiver LOS must fulfil the relations: According to this geometry, the free electron density distribution can be described with a two-dimensional function, N (r, ϕ), so that Eq. ( 1) gives rise to: ) where dr /cos z = dγ , and h 1 ∼ 100 km and h 2 ∼ 1,000 km are the nominal lower and upper boundaries of the iono-sphere.At this point, it is worth stressing the fact that not subindexes are added to vT EC or sT EC notation in Eq. ( 6) in order to indicate that these are exact values computed from the actual electron density distribution.Using Eqs. ( 2) and ( 6), the following expression can be obtained for the error in the conversion from sT EC to vT EC L (hereafter, the conversion error, ε), when an ILH h L is used, and for a given LOS that passes through an IPP with geographic latitude ϕ L , observed from a GNSS receiver at geographic latitude ϕ R : The mean and the standard deviation (Kalton 1983) of the conversion error for all the observing LOS provide, respectively, estimates of the systematic and non-systematic components of the conversion error for a given ILH and geographic latitude: Just to provide an example, Fig. 3 shows the variation of the conversion error with the IPP latitude for different ILH, for a GNSS receiver at −20 • of geographic latitude, for high solar activity and March.The values of the mean and the standard deviation are also presented in the legend of this figure.This example serves to illustrate the sensitivity of the conversion error to the ILH.In particular, it can be appreciated that there exists an ILH (between 400 and 450 km) that reduces to zero the mean of the conversion error, but there does not exists a height that reduces to zero its standard deviation.Now, a chain of seven GNSS receivers along the −60 • meridian was simulated for this assessment.The observing Figure 4 shows these estimates for high and low solar activity, and for March and June.These few cases show that the mean is zero for an ILH, h L ,0 , that varies with the solar activity and month.
Figure 5 presents a complete description of the variability of the "optimal IHL", h L ,0 , i.e., the one that fulfills the condition μ ε (h L ) = 0, and σ ε h L ,0 with solar activity and month.Both estimates show similar variability patterns: (i) the values increase with the solar activity; (ii) the dependence on the solar activity is greater in December than in June solstice; and (iii) the greatest values take place in both equinoxes, the lowest values take place in the December solstice, and intermediate values take place in the June solstice.The patterns observed in h L ,0 are consistent with the variability that presents the hmF2 and the foF2 and has been registered by ionosondes (Rishbeth 2000;Li and Yu 2003) that are known as the annual and semi-annual anomalies of the ionosphere (Mendillo et al. 2005;Azpilicueta and Brunini 2011).
As can be seen in the upper panel of Fig. 5, the range of the h L ,0 values for the December solstice are approximately 340 and 410 km; for the equinoxes are approximately 390 and 410 km; and for the June solstice are approximately 370 and 390 km.The standard deviation, σ ε h L ,0 , varies between approximately ±0.5 and ±2.5 TECu during the December solstice; ±1.1 and ±2.5 TECu during the equinoxes; and ±1.0 and ±1.8 TECu during the June solstice.

Assessment of the vTEC and IFB errors
The conversion error discussed in the previous section propagates to the unknowns of Eq. ( 3) when their values are estimated from the GNSS observations and causes an incorrect estimation of both, the x 0 , . .., x n parameters of the vT EC L function and the IFB, β.In order to evaluate the error propagation, the following function is proposed to represent the latitudinal variation of the vT EC along the −60 • meridian, from −50 • to +30 • of geographic latitude, at 14 LT: where P l,0 (sin ϕ) are the fully normalized associated Legendre functions of degree l and order zero, and x i , i = 0, . . ., n, is a set of constant expansion coefficients to be determined.The use of Legendre functions is based on the fact that the common choice in global TEC mapping as function of latitude, longitude and time is the use of surface harmonics, based on Legendre functions that are orthogonal over a sphere (Schaer 1999).Since the ionospheric scenarios The first experiment performed in this section consists of determining the expansion coefficients in order to fit, with the Least Squares criterion, the vT EC computed with the NeQuick model (using Eq. 6a).This experiment confirms that an expansion of maximum degree L = 12 provides an excellent representation of the NeQuick vTEC, with a standard deviation for the residuals of the LSQ adjustment well below ±0.1 TECu for march and high solar activity, the worst case scenario, with lower standard deviation values for the other cases.
Once verified that the function proposed in Eq. ( 9) with maximum degree L = 12 is well suited for representing the NeQuick vT EC, a second experiment is performed in order to fit the same function to the sT EC computed with the NeQuick model: where η is added to account for the MF errors.The standard deviations, σ η (h L ), obtained after the Least Squares fit are: ±5.3 TECu for h L = 350 km, ±3.2 TECu for h L = 400 km, ±3.7 TECu for h L = 450 km, and ±5.7 TECu for h L = 500 km (Table 1).These larger standard deviations (with respect to the ones obtained in the previous Table 1 Standard deviations obtained after the LSQ adjustment for Eq.(10a) that does not include the IFB, σ (1) η (h L ), and for Eq. ( 11) that includes the IFB, σ experiment) provide the first assessment regarding the effects of the conversion error.
The curves of different colours in the left-hand side panel of Fig. 6 show the vT EC L computed after fitting the expansion coefficients x 0 , . . ., x n for different ILH; while the black curve shows the "true" vT EC computed with the NeQuick model.
In order to characterize the errors caused by the MF in the vT EC L estimation, the differences between this estimate and the "true" vT EC is computed: This error can be characterized by its mean and standard deviation: These estimates are shown in the right-hand side panel of Fig. 6 for high solar activity and March.According to this plot, the values of μ vT EC (h L ) and σ vT EC (h L ) range from -3.6 ±1.4 to +3.1 ±1.6 TECu, reaching their minimum for an ILH of 425 km. Figure 7 presents a complete description of the variability with solar activity and month of the standard deviation, σ vT EC h L ,0 , that corresponds to the ILH h L ,0 that fulfils the condition μ ε (h L ) = 0.The values in that figure ranges from ±0.5 to ±2.0 TECu.
The last experiment performed in this section consists in adding the IFB, β, to Eq. (10a) and to estimate their values together with the expansion coefficients x 0 , . . ., x n : The inclusion of additional unknowns in the Least Squares adjustment reduces (with respect to the ones obtained in the Since the sT EC computed with the NeQuick model is not affected by any IFB, the "true" value of the IFB is β = 0; in other words: any deviation from zero in the estimated β unknowns must be interpreted as an error, β, caused by the propagation of the conversion errors to the estimated unknowns.The lines of different colours in the left panel of Fig. 9 show the β values obtained for different ILH, for high solar activity and March, after solving Eq. ( 11).
As the other estimates previously discussed in this section, the β error is characterized by its mean, μ β (h L ), and standard deviation, σ β (h L ), for the seven GNSS receivers that were simulated in this study.These estimates are shown in the right panel of Fig. 9 for high solar activity and March.
It should be noted that the ILH, h L ,0 , that cancels the mean vT EC L error is, in general, different from the ILH, h L ,0 , that cancels the mean IFB error.The difference h L ,0 = h L ,0 − h L ,0 is represented in the upper panel of Fig. 10 for different solar activity conditions and months.From March to September h L ,0 varies from 5 to 12 km, almost independently of the solar activity.From September to March h L ,0 is strongly dependent on the solar activity: for low solar activity varies from −1 to −12 km, while for high solar activity varies from +3 to +4 km.The bottom panel of Fig. 10 shows the standard deviation of the IFB error σ β h L ,0 , for different solar activity conditions and months.The values range from approximately ±0.7 to ±3.5 TECu, the greatest contribution to this variability coming from the southern GNSS receiver located at -40 o of geographical latitude.σ β h L ,0 is greater for high than low solar activity, and greater for equinoxes than solstices.

Conclusions
According to the numerical simulation performed in this work, there does not exist a unique ILH that reduces to zero the conversion error (i.e., conversion from slant to vertical TEC) for the whole region study in this research (i.e., the Central and the northern part of the South American continents).For a given solar activity and season, there exists an ILH that reduces to zero the mean conversion error, but still remains a latitude-dependent conversion error.
The ILH that reduces to zero, the mean conversion error varies with the solar activity between approximately 340 and 410 km during the December solstice, 390 and 410 km during the equinoxes, and 370 and 390 km during the June solstice.When that ILH is used, the latitude-dependent conversion error reaches a standard deviation that increases with the solar activity and varies approximately between ±0.5 and ±2.5 TECu during the December solstice; ±1.1 and ±2.5 TECu during the equinoxes; and ±1.0 and ±1.8 TECu during the June solstice.
The latitude-dependent conversion error propagates to the parameters of the model used to represent the latitudinal  Based on this study, the use of a fixed 'arbitrary' IHL regardless of the season and solar activity is not correct and it should be choosen considering the annual and semiannual variability that present with month and solar condition.Figure 5 shows the ILH that reduced the mean errors to zero, at least within the conditions simulated for this study, and the observed variability patterns could help in the choice of the ILH.
Besides, the ILH that cancels the mean vertical TEC error is different from the one that cancels the mean IFB error.From March to September, the difference between both heights varies from 5 to 12 km, almost independently of the solar activity.From September to March, that difference is strongly dependent on the solar activity: for low solar activity varies from -1 to -12 km, while for high solar activity varies from +3 to +4 km.If the height that cancels the mean vertical TEC error is used, the IFB error reaches a standard deviation that varies approximately from ±0.7 to ±3.5 TECu, depending on the solar activity and month.It is greater for high than low solar activity, and greater for equinoxes than solstices.

Fig. 1 Fig. 3
Fig. 1 cuts of the electron density distribution computed using the Nequick model along the −60 • meridian, from 100 to 1,000 km of height above the Earth's surface (the solid line at 6,720 km

Fig. 5
Fig.5Upper panel shows the variation of the "optimal IHL", h L ,0 , the height layer that fulfills the condition μ ε (h L ) = 0, in km, as a function of month of the year (x-axis) and solar activity (y-axis).Bottom panel shows the variation of σ ε h L ,0 in TECU, as a function of month of the year (x-axis) and solar activity (y-axis)

Fig. 6 Fig. 7
Fig. 6 vT EC L estimated from Eqs. (10) for ILH of 350, 400, 450 and 500 km (different colours) and "true" vT EC (black; left panel, in TECU); mean and standard deviation of the differences between vT EC L and vT EC(right panel, in TECu); for high solar activity and March

Fig. 8
Fig. 8 vT EC L estimated from Eq. (11) for ILH of 350, 400, 450 and 500 km (different colours) and "true" vT EC (black; left panel, in TECU); mean and standard deviation of the differences between vT EC L and vT EC(right panel, in TECu); for high solar activity and March