An SBAS Integrity Model to Overbound Residuals of Higher-Order Ionospheric Effects in the Ionosphere-Free Linear Combination

The next generation of satellite-based augmentation systems (SBAS) will support aviation receivers that take advantage of the ionosphere-free dual-frequency combination. By combining signals of the L1 and L5 bands, about 99% of the ionospheric refraction effects on the GNSS (Global Navigation Satellite Systems) signals can be removed in the user receivers without additional SBAS corrections. Nevertheless, even if most of the negative impacts on GNSS signals are removed by the ionospheric-free combination, some residuals remain and have to be taken into account by overbounding models in the integrity computation conducted by safety-of-live (SoL) receivers in airplanes. Such models have to overbound residuals as well, which result from the most rare extreme ionospheric events, e.g., such as the famous “Halloween Storm”, and should thus include the tails of the error distribution. Their application shall lead to safe error bounds on the user position and allow the computation of protection levels for the horizontal and vertical position errors. Here, we propose and justify such an overbounding model for residual ionospheric delays that remain after the application of the ionospheric-free linear combination. The model takes into account secondand third-order ionospheric refraction effects, excess path due to ray bending, and increased ionospheric total electron content (TEC) along the signal path due to ray bending.


Introduction
Airplanes need to rely on ground-based radio systems such as the instrument landing systems (ILS) or satellite-based augmentation systems (SBAS) to be able to perform precise landing operations independent of the weather conditions. Given the fact that SBAS, similar to the European Geostationary Overlay System (EGNOS) or the American Wide Area Augmentation System (WAAS, allow an efficient use of airspace without maintenance or calibration cost for airports, more and more airports are already taking advantage of this technology based on GNSS (Global Navigation Satellite System). (For instance, according to the official website of the FAA (Federal Aviation Administration, already 4053 WAAS procedures are currently (21 May 2020) implemented in 1956 US airports.).
However, ionospheric effects on trans-ionospheric GNSS signals are significant and should not be ignored. The ionospheric propagation effects degrade the true range measurements, whereas the ionospheric irregularities may cause signal loss due to rapid fluctuation of the signals, which is also known as ionospheric scintillation. Signal fluctuations may cause cycle slips; however, if these cycle slips are detected, precise navigation can be achieved under scintillation conditions [1]. Tiwari and Strangeways [2] presented an approach to mitigate the effect of ionospheric scintillation on GNSS users in the European region using total electron content at 1 Hz rate. Demyanov et al. [3] found degradation in the single-frequency WAAS services during geomagnetic storm conditions. Fortunately, the next generation will not only support single-frequency receivers but also aviation receivers that benefit from the ionospheric-free dual-frequency combination. By combining signals of the L1 and L5 bands, future users will be able to remove about 99% of the ionospheric refraction effects on the GNSS signals instantaneously in the receivers on board the airplane. Thus, this technique will help to significantly increase the service performances of SBAS, in particular for polar and equatorial latitudes where geomagnetic storms or strong ionospheric gradients accompanied by phase and amplitude scintillations may lead to partial losses of the service availability. However, even if most of the negative impacts on GNSS signals are removed by this combination, some residuals remain. They will need to be taken into account in the protection level computation conducted within aircraft receivers by the use of so-called "overbounding models". Such models overbound residuals also, which result from most rare extreme ionospheric events such as the famous "Halloween Storm" and should thus include the tails of the error distribution. Their application shall lead to safe error bounds on the user position and allow the computation of protection levels for the horizontal and vertical position errors.
Brunner and Gu [4], Bassiri and Hajj [5], and Jakowski et al. [6] were among the first to compute higher-order ionospheric effects in GPS dual frequency combinations. Later, Kedar et al. [7], Fritsche et al. [8], and Hernández-Pajares et al. [9] studied the second-order effect on receiver positions, satellite clocks, and orbits. They found a global southward displacement of several millimeters in the station positions. Hoque and Jakowski [10][11][12][13] derived different approximation formulas to correct the second-and third-order terms in the refractive index, errors due to geometric bending, and differential slant total electron content (TEC) bending effects on geometric range estimation. A comprehensive summary of higher-order effects on geodetic quantities considering all relevant higher-order contributions can be found in Hernández-Pajares et al. [14].
Hadas et al. [15] investigated the impact of the second-and third-order terms, geometric bending effects on satellite orbits, satellite clock corrections, receiver positions, etc. They found that the satellite-related products captured most of the impact of higher-order corrections with the magnitude up to 2 cm for clock corrections. Zhang et al. [16] investigated the impact of higher-order ionospheric effects on tropospheric delay estimation using GNSS observations from about 100 worldwide ground stations. They found that the tropospheric zenith wet delay (ZWD) estimation is affected up to 20 mm during ionospheric and geomagnetic storms. Cai et al. [17] used a quad-constellation (GPS, GLONASS, BDS and Galileo) precise point positioning (PPP) approach for investigating higher-order ionospheric impact in the position domain. They found that higher-order ionospheric terms can affect the 3D position solutions of the quad-constellation PPP up to 6 mm.
Here, we propose and justify a simple overbounding model that bounds ionospheric-free residual delays and can be used for the calculation of the protection levels of an airplane. The model is particularly designed for the "Minimum Operational Performance Standards for Galileo-Global Positioning System-Satellite-Based Augmentation System Airborne Equipment" [18]. The model takes into account the second-and third-order ionospheric refraction effects, excess path due to ray bending, and increased ionospheric TEC along the signal path due to ray bending (the impact of phase and code noises are not part of the model).
The design of the model is very much driven by constrains of the aviation users. The requirements are: • the model shall basically overbound any ionospheric-free higher-order residual everywhere for static and dynamic users; • the model shall be easy to implement in an aviation receiver and self-standing in terms of not relying on external parameters; • the model does not have to provide realistic data in terms of the real characteristics of the ionosphere, as long as it is overbounding the residuals and the maximum values of the model do not jeopardize the error budget of the future aviation user. The overbounding model is developed by a method that can be seen as an educated guess, combining theoretical knowledge with a number of extreme assumptions on some physical parameters. This approach is chosen because the model will be hard coded in the user receiver, and we have to assure that our model overbounds anomalous ionospheric delays even for the future. This is quite different to single-frequency SBAS. There, the SBAS itself generates the ionospheric error bounds for the user with respect to the broadcast user corrections. This approach allows the post-adoption of the models for the ionospheric threat space based on continuous monitoring activities. For a more detailed overview on SBAS safety assessment guidance related to anomalous ionospheric conditions, we refer to [19].
Moreover, one has to keep in mind that the alert limits for positioning errors for CAT-I landing operations are 40 m for the horizontal errors and 12 m for the vertical errors. Given the fact that the first-order ionospheric effects, which are substantially contributing to the error budget in the singlefrequency case are removed, statistical uncertainty can be compensated by putting additional safety margins on the assumptions that are the foundation of the overbounding model.
The baseline criteria to design an overbounding model according to the above listed requirements are: • to design the model, using an almost deterministic approach grounded on physical equations; • to define input parameters for the equations that are beyond what has been reported in the literature or close to the most extreme values (to estimate representative values, with a very low occurrence probability all over the globe and under all solar conditions would be most likely impossible for an environment as complex as the ionosphere); • to put additional margins on the input parameter, taking into account that there may be rare events that are never observed by scientists (such as what has been named a Black Swan by the risk analyst Nassim Nicolas Taleb [20]), but ensuring that the conservatism of the model does not destroy the overall bounding; • to validate the model in a multi-dimensional approach to ensure that the mapping function and the approximations, inherent in some of the equation used, are not destroying the overbounding of the model or leading to modeling errors.
The first section gives a brief overview of the impact of ionospheric refraction effects on GNSS signals, followed by a mathematical description of the residuals in the ionospheric-free combination. After a short section on the impact of carrier smoothing on the residuals, we detail the design and justification of the overbounding model according to the strategy that has been detailed above. The subsequent sections describe the ray-tracing results, verifying that the used approximation to compute the third-order residuals and imprecisions introduced by the use of the mapping function do not lead to a violation of the conservativeness of the model. Finally, the paper is closed by a brief summary.

Ionospheric Refraction Impacts on Signal Propagation
Ionospheric effects on radio signals can be mathematically formulated as the expansion of the Appleton-Hartree formula [21]. In consideration of different propagation velocities for a single phase carrier and a wave group, the observables for GNSS code (Ψ) and carrier-phase (Φ) observations can be expressed in the dimension of length [22] by the following equations: Remote Sens. 2020, 12, 2467 4 of 17 Here, the index i refer to the GNSS carrier frequency, e.g., L1 or L5. The parameters Ψ and Φ are the measured code and carrier phase ranges, respectively. λ is the corresponding wavelength, and n is the integer ambiguity. ρ corresponds to the geometrical distance plus clock errors and non-ionospheric delays. The quantity d b(length) i denotes the error contribution of the excess path length due to geometric bending of the signal within the ionosphere. In both equations, the first-order term is related to as f 2 , the second-order term is related to as f 3 , and the third-order term is related to as f 4 . Other error sources such as carrier phase ambiguities, satellite and receiver instrumental biases, multipath effects, and measurements errors are omitted for simplicity in the formulas. The ionosphere-related parameters p, q, and t in (1) and (2) are given by: t = 2437 n 2 e ds + 4.74 × 10 22 n e B 2 1 + cos 2 Θ ds.
Here, n e denotes the electron density and s denotes the signal propagation path. TEC LoS is the total electron content (TEC) along the line of sight between receiver and satellite. TEC bend is the TEC deviation due to the difference between the real curved ray path and the ray path along the straight line of sight. TEC bend is not removed by the ionospheric-free combination, because the ray path is individual for every frequency. B is the magnetic induction, and Θ is the angle of intersection between the Earth's magnetic field vector and the propagation direction of the transmitted signal. The International System of Units (SI) is used when computing each variable.
To determine the geometric bending term b(length) and TEC bend , we use the formulas proposed by Hoque and Jakowski [13]: In (6), b1 = 2.9344×107 and b2 = 0.8260 are model coefficients, and β is the elevation angle in radians. Here, the signal frequency f is in MHz and b(length) is in centimeters. Equation (7) computes TEC bend in TECU, H is the atmospheric scale height in km, hmF2 is the maximum ionization height in km, and TEC LoS has to be in TECU, β again is in radians, and f isin GHz.
According to Bassiri and Hajj [5], the second term of the ionospheric parameter t leads just to a sub-millimeter range bias for GHz frequencies, and the same appears in terms of higher order. Consequently, these terms will be neglected in the following sections [8,23,24].

Residuals in the Ionospheric-Free Combination
The ionospheric-free combination for the code and carrier phase (in units of length) on f i and f j frequency (f i > f j ) is given by: Remote Sens. 2020, 12, 2467 5 of 17 where Here, ∆s b(TEC) denotes residual excess TEC caused by the difference in ray bending of the two combined signals indexed by i and j. ∆s 2 and ∆s 3 are the dual frequency second-and third-order residuals, ∆s b(length) is the residual excess ray-path length caused by the difference in ray bending of the two signals, and A ifc is the ionospheric-free linear combination of carrier phase ambiguities.

Impact of Carrier-Smoothing
Datta-Barua et al. [24] have shown that the higher-order group delay residuals, and to a much lesser extent the higher-order phase delay residuals, contribute to the smoothing error in a carrier smoother (Hatch filter) designed according to RTCA MOPS (Radio Technical Commission for Aeronautics Minimum Operational Performance Requirements) [25]. Due to the mathematical formulation of the carrier smoother, only changes in the higher-order phase delay residuals appear as phase error rates in the filtering process.
The authors have also shown that for a 100 s Hatch-filter time constant, a 20 cm higher-order group delay will exhibit a 20 cm smoothing error. This impact can be reduced to a few centimeters with a longer filter time constant, as they have shown for a 2-h filter. Thus, we consider the impact of the filter time for the higher-order group delays as almost negligible for an aviation receiver with filter times between 100 and 600 s.
With respect to the higher-order phase delay residuals that appear as phase error rates in the carrier smoother, the authors have shown that the higher-order phase error rate has no significant effects in either the 100 s or two-hour filter estimates. Thus, we consider the impact of the filter time for the higher-order phase delays as negligible for an aviation receiver.

Design and Justification of the Overbounding Model
The SBAS integrity equation has been introduced by Walter et al. [26]. It describes the position error distribution after the application of the SBAS corrections. A protection level is a statistical bound of the horizontal or vertical position error. One input for the computation of protection levels is the standard deviation of a zero mean normal distribution of the ionospheric residual error after the ionospheric-free combination (σ DF_residual ). The σ DF_residual has to be designed in a way that the underlying normal distribution overbounds any error up to 10 −7 . As mentioned in the introduction, the aim is to define an overbounding model that is implementable in an aviation receiver. The latter leads to the following requirements: 1.
Only parameters that are available in an aviation receiver shall be used; 2.
Simple mathematical formulation to ease the implementation in a receiver; 3. Flexibility (a model update is just a change of a few configuration parameters).
Remote Sens. 2020, 12, 2467 6 of 17 Taking into account these requirements, the mathematical formulation of the residuals and the conclusions on the impact of the Hatch filter implemented in an aviation receiver, we have taken some design choices and assumptions for the modeling approach: • The final overbounding model shall be elevation-dependent; • The main driver for the model design is the TEC; • The dominating factor for the higher-order residuals are the group delay residuals; due to their low impact (see section above), potential phase delay residuals are not taken into account; • The filter constant of the Hatch filter used in the data pre-processing of an aircraft is in the order of 100 s. Therefore, the filter is assumed to have no limiting impact on the residuals.

Calculation of the Overbounding Model
As shown by (4)-(6), the physical parameters driving the magnitude of the ionospheric-free residuals are the total electron content TEC, the electron density ne, and the magnetic induction B. Our strategy to design an overbounding model is a five-step approach: 1.
Take worst-case assumptions on the geometry and magnetic induction B; 2.
Use approximations that link the electron density in (5) to TEC; 3.
Define and adjust the overbounding model function and their parameter according to the data; 5.
Perform a ray-tracing, to ensure that our approximations and the use of the ionospheric mapping function are not violating the conservatism of the approach.

Worst-Case Magnetic Induction
As apparent from (4), the geomagnetic induction is one of the impact parameters for second-order ionospheric effects. The International Geomagnetic Reference Field (IGRF) model 11th generation [28] is used to compute gridded maps of the geomagnetic induction at the Earth's surface and at an ionospheric height of 350 km with 2.5 degree latitude and longitude resolution for an arbitrary date 1 January 2020. Analyzing these values, we found that the maximum geomagnetic induction is about 0.67 and 0.56 Gauss at the Earth's surface and at an ionospheric height of 350 km, respectively. As a worst-case value for the geomagnetic induction, we choose 0.6 G, which reflects the strength of the geomagnetic field at ionospheric altitudes in the Polar Latitudes. Given that the magnetic induction is maximal in the Polar Regions, this value is seen as a representative upper bound for the global geomagnetic field strength.
To compute the overbounding model we make the following further assumptions: • The magnetic field is always constant (worst case) along the ray-path; • The azimuth is always set in a way that the magnetic induction is maximal (signal path parallel to magnetic field lines)

Linking Electron Density to TEC
As a first step to link the electron density to the TEC, we used the approximation of Brunner and Gu [4]: Here, the integral of the squared electron density along the line of sight is replaced by the peak electron density Nmax, the total electron content along the line of sight (TECLoS), and a shape parameter η. The latter is set to 0.66 as an appropriate value to account for different electron density distributions [8,22]. Finally, to find an appropriate value for N max , we applied the model introduced by Pireaux et al., in 2010 [29]: The model links Nmax to the vertical TEC (VTEC). To complete the model to overbound the ionospheric-free residuals, the remaining task is to find extreme cases of TEC that fit the assumptions of the integrity risk.

Selection of Extreme TEC Values
The first step to identify extreme TEC values is the analysis of the 42 most severe geomagnetic storm events that occurred since 1994 (selected by the geomagnetic indices Ap and kp), to identify dates and regions where the highest TEC values occurred. The TEC values and their associated locations are derived from post-processed daily 2-hourly CODE (Center for Orbit Determination Europe) TEC maps. The spatial resolution of these maps is 2.5 • and 5 • along latitude and longitude, respectively. A part of the outcome of this analysis is shown in Table 1. It lists the four dates and corresponding locations where the highest TEC values are found.  For these stations, we processed the TEC from the GPS observations for the complete days with a time resolution of 30 s. An example of the output of this process is depicted in Figure 1 Table 2.  The maximum TEC values that we found for the 5 stations for 29 October 2003 and 30 October 2003 through this procedure are listed in Table 2. Given the fact that the highest TEC values occurred on 30 October 2003 at the stations LHUE, HNLC, and SIO3, we further limited our analysis of the data of these three stations on only this date.
To derive a TEC value that represents a probability beyond 10 −7 from this data, we looked for an appropriate statistical prediction model. First, we crosschecked the data against the following probability distributions: Normal, Exponential, Extreme value, Lognormal, and Weibull. Examples of the outcomes of this analysis for the Normal and Extreme value distributions are displayed in Figures 2-4. We created a normal probability plot (e.g., Figure 2a) comparing the distribution of the data to the normal distribution. The broken reference line indicates whether the data follows a normal distribution or not. If the data points appear along the reference line, it confirms that the sample data has a normal distribution (see Figure 3a). If the distribution is not normal, it will introduce curvature in the data plot (see Figure 2a). Similarly, we created an extreme value probability plot (e.g., Figure 2b) comparing the distribution of the data to the extreme value distribution, and the broken reference line indicates whether the data follows an extreme value distribution or not. We used MATLAB built-in functions for obtaining probability distribution [30].
In a second step, we calculated the TEC from the inverse of the cumulative density function (CDF) of our model distributions. The corresponding cumulative distribution functions are given in Figure 5, and the final results of this procedure are depicted in Table 3. The latter lists the maximum TEC values obtained from the inverse CDF for a probability bound of 10 −7 .
Remote Sens. 2020, 12, 2467 8 of 17 Given the fact that the highest TEC values occurred on 30 October 2003 at the stations LHUE, HNLC, and SIO3, we further limited our analysis of the data of these three stations on only this date.
To derive a TEC value that represents a probability beyond 10 −7 from this data, we looked for an appropriate statistical prediction model. First, we crosschecked the data against the following probability distributions: Normal, Exponential, Extreme value, Lognormal, and Weibull. Examples of the outcomes of this analysis for the Normal and Extreme value distributions are displayed in Figures 2-4. We created a normal probability plot (e.g., Figure 2a) comparing the distribution of the data to the normal distribution. The broken reference line indicates whether the data follows a normal distribution or not. If the data points appear along the reference line, it confirms that the sample data has a normal distribution (see Figure 3a). If the distribution is not normal, it will introduce curvature in the data plot (see Figure 2a). Similarly, we created an extreme value probability plot (e.g., Figure 2b) comparing the distribution of the data to the extreme value distribution, and the broken reference line indicates whether the data follows an extreme value distribution or not. We used MATLAB built-in functions for obtaining probability distribution [30].     In a second step, we calculated the TEC from the inverse of the cumulative density function (CDF) of our model distributions. The corresponding cumulative distribution functions are given in Figure 5, and the final results of this procedure are depicted in Table 3. The latter lists the maximum TEC values obtained from the inverse CDF for a probability bound of 10 −7 .    Only after comparing the maximum TEC values observed (first column of Table 3) to the extrapolated TEC values for a probability bound of 10 −7 obtained from inverse CDF analysis, we found that the Normal and Extreme Value distributions show acceptable results, whereas the Exponential, Lognormal, and Weibull distributions give TEC values far beyond any physical explanation. Since the probability distribution tests of the station SIO3 show a bad fitting performance of the models for the Normal and Extreme Value distributions (Figure 2, top and panels), we took the data from SIO3 out of our further considerations. The tests for the stations LHUE and HNLC show a reasonable fitting performance for both Normal and Extreme Value Distribution, as shown in Figures 3 and 4. To conclude the selection process, we decided to take the higher of the two remaining values. Therefore, 354.0 TECU obtained from HNLC is the selected reference, and we determined 360 TECU as the maximum TEC to design our overbounding model. Obviously, this estimation has not led to a precise TEC estimate, but to an overestimated TEC as required for the overbounding model. In a second step, we calculated the TEC from the inverse of the cumulative density function (CDF) of our model distributions. The corresponding cumulative distribution functions are given in Figure 5, and the final results of this procedure are depicted in Table 3. The latter lists the maximum TEC values obtained from the inverse CDF for a probability bound of 10 −7 .  Only after comparing the maximum TEC values observed (first column of Table 3) to the extrapolated TEC values for a probability bound of 10 −7 obtained from inverse CDF analysis, we found that the Normal and Extreme Value distributions show acceptable results, whereas the Exponential, Lognormal, and Weibull distributions give TEC values far beyond any physical

The Overbounding Model
To compute the overbounding model, one has to sum all residuals of the ionospheric-free combination for the worst-case conditions as described in the previous sections. In our approach, we use (8)-(10) with a maximum TEC value of 360 TECU, the mapping function as described in RTCA/DO-229D, a constant magnetic induction of 0.6 G, which is obtained at the magnetic poles, and an elevation range from 2 to 90 degrees. The choice of 0.6 G for all latitudes has the advantage that it will put additional margin (in terms of increasing the overbound) to lower latitudes, where the magnetic induction is lower. The results of the computations are depicted in Figure 6. It shows the elevation-dependent curves for the sum of all residuals, the sum of second-and third-order effects, the excess path, excess TEC, and the final overbounding model. To conclude the selection process, we decided to take the higher of the two remaining values. Therefore, 354.0 TECU obtained from HNLC is the selected reference, and we determined 360 TECU as the maximum TEC to design our overbounding model. Obviously, this estimation has not led to a precise TEC estimate, but to an overestimated TEC as required for the overbounding model.

The Overbounding Model
To compute the overbounding model, one has to sum all residuals of the ionospheric-free combination for the worst-case conditions as described in the previous sections. In our approach, we use (8), (9), and (10) with a maximum TEC value of 360 TECU, the mapping function as described in RTCA/DO-229D, a constant magnetic induction of 0.6 G, which is obtained at the magnetic poles, and an elevation range from 2 to 90 degrees. The choice of 0.6 G for all latitudes has the advantage that it will put additional margin (in terms of increasing the overbound) to lower latitudes, where the magnetic induction is lower. The results of the computations are depicted in Figure 6. It shows the elevation-dependent curves for the sum of all residuals, the sum of second-and third-order effects, the excess path, excess TEC, and the final overbounding model. The final model (see Figure 7) is obtained by dividing the sum of all residuals by a factor of 5.33 to convert the standard deviation for the 10 −7 probability to a one σ standard deviation and to fit a selected function that fulfills the requirement to be easy implementable in GNSS receiver software. This approach leads to the following overbounding model for the variance of the ionospheric-free residuals.
"When applying the ionospheric-free dual-frequency L1/L5 combination, the residual error (including ionospheric higher-order effects, ray bending, and excess TEC) in units of meters is The final model (see Figure 7) is obtained by dividing the sum of all residuals by a factor of 5.33 to convert the standard deviation for the 10 −7 probability to a one σ standard deviation and to fit a selected function that fulfills the requirement to be easy implementable in GNSS receiver software. This approach leads to the following overbounding model for the variance of the ionospheric-free residuals.
"When applying the ionospheric-free dual-frequency L1/L5 combination, the residual error (including ionospheric higher-order effects, ray bending, and excess TEC) in units of meters is described by the normal distribution N(0, σ 2 i,UIRE ), where: and El i is elevation angle in degrees for the i th satellite. The model is valid down to an elevation angle of 3 degrees." Remote Sens. 2020, 12, 2467 11 of 17 Figure 7. The final overbounding model.

Model Justification by Ray-Tracing
As a final step for the justification of the model, we applied a ray-tracing technique. This method validates that the used approximations in computing different higher-order residuals and imprecisions introduced by the use of the mapping function do not lead to a violation of the conservativeness of the model. The approach is based on a two-dimensional ray-tracing technique developed by Hoque and Jakowski [12]. The anisotropy of the ionosphere, i.e., the effect of the Earth's magnetic field on the signal propagation, has been taken into account. For this, the International Geomagnetic Reference Field (IGRF, [28]) model is applied when calculating geomagnetic field vectors at numerous points along incoming ray paths. The interactions among the signal, ionospheric electron density, and geomagnetic field lines are approximated by the Appleton-Hartree formula [31]. The ionosphere is assumed to be composed of numerous thin spherical layers in each of which the medium is homogeneous.

Modeling of the Earth's Ionosphere and Plasmasphere
The electron density distribution of the Earth's ionosphere is modeled by adding multiple Chapman layers superposed by an exponential decay function representing the plasmasphere. The Chapman layer [32] describes the electron density Ne distribution as a function of height (h) by is the ionospheric electron density and (ℎ) is the plasmaspheric electron density, χ is the sun's zenith angle, N0 is the peak electron density observed at the altitude h0 at zenith (χ = 0), Figure 7. The final overbounding model.

Model Justification by Ray-Tracing
As a final step for the justification of the model, we applied a ray-tracing technique. This method validates that the used approximations in computing different higher-order residuals and imprecisions introduced by the use of the mapping function do not lead to a violation of the conservativeness of the model. The approach is based on a two-dimensional ray-tracing technique developed by Hoque and Jakowski [12]. The anisotropy of the ionosphere, i.e., the effect of the Earth's magnetic field on the signal propagation, has been taken into account. For this, the International Geomagnetic Reference Field (IGRF, [28]) model is applied when calculating geomagnetic field vectors at numerous points along incoming ray paths. The interactions among the signal, ionospheric electron density, and geomagnetic field lines are approximated by the Appleton-Hartree formula [31]. The ionosphere is assumed to be composed of numerous thin spherical layers in each of which the medium is homogeneous.

Modeling of the Earth's Ionosphere and Plasmasphere
The electron density distribution of the Earth's ionosphere is modeled by adding multiple Chapman layers superposed by an exponential decay function representing the plasmasphere. The Chapman layer [32] describes the electron density N e distribution as a function of height (h) by N I e (h) = N 0 exp(0.5(1 − z − secχ · exp(−z))).
is the ionospheric electron density and N P e (h) is the plasmaspheric electron density, χ is the sun's zenith angle, N 0 is the peak electron density observed at the altitude h 0 at zenith (χ = 0), and H is the atmospheric scale height. It can be shown that the maximum ionization Nm = N 0 √ cosχ occurs at the height hm = h 0 + Hln(secχ) [21]. The plasmaspheric density N P e (h) is modeled by an exponentially decreasing function of height where n p is the plasmaspheric basic density of electrons and H P (10,000 km) is the mean scale height of the plasma density [33]. The plasmaspheric density contribution below the peak density height hm is set to zero.
The atmospheric scale height H, peak density Nm, and peak height hm practically define the shape of the Chapman layer; thus, by tuning these parameters, we can generate numerous electron density profiles with varying shapes. The expression for vertical TEC can be derived by integrating Equation (14) with respect to the height and written as [11] In order to simulate a broader variety of profile shapes, we superposed three Chapman layers representing the ionospheric E, F1, and F2 layers. The values of shape parameters H, Nm, and hm used for individual layers are determined by extensive literature research (e.g., [21,[34][35][36][37][38][39]). We found that they could vary largely depending on geophysical conditions, which again depend on geographic/geomagnetic location as well as diurnal, seasonal, and solar cycle variation. Based on literature study, we set separate ranges for peak electron density heights hmF2, hmF1, hmE, and scale heights HF2, HF1, and HE of ionospheric F2, F1, and E layers as follows: 250 ≤ hmF2 ≤ 550 km, 60 ≤ HF2 ≤ 90 km 180 ≤ hmF1 ≤ 220 km, 40 ≤ HF1 ≤ 60 km 100 ≤ hmE ≤ 120 km, 5 ≤ HE ≤ 15 km Knowing the scale height and TEC contribution of the F2, F1, and E layers, the corresponding peak electron densities NmF2, NmF1, and NmE can be determined using Equation (16) by where VTECF2, VTECF1, and VTECE are the VTEC contributions from individual layers. The absolute maximum VTEC is determined as 360 TECU in a previous section. In order to describe a broader variety of geophysical conditions, we assume that about 5-20% (randomly fixed) of the ionospheric electron contribution comes from the F1 layer, about 0-5% (randomly fixed) comes from the E layer, and the rest of the electron contribution comes from the F2 layer (for more details, see Kelley [35] and Norman [37]). The solar zenith angle χ is varied within the range 0-80 • when modeling each layer, and values above 80 • are ignored to avoid complexity in the ionospheric model (Davies [36]). We further assumed that the plasmaspheric contribution could vary from about 5% to 10% (randomly fixed) of the total contribution.

Justification of Mapping Function
The ray-tracing results are used to verify the accuracy of the ionospheric mapping function. The ionospheric mapping function (MF) is used for vertical to slant TEC conversion and vice versa. There are several standards for the mapping function depending on the application (Figure 8). Here, we compare the MF recommended in MOPS 229D (see 419, 2008), GPS ICD (see GPS Joint Program Office 2004), and the modified single-layer model mapping (MSLM) given by Dach et al. [40] with ray-tracing results. The MF algorithms can be found in the given references.

Justification of Mapping Function
The ray-tracing results are used to verify the accuracy of the ionospheric mapping function. The ionospheric mapping function (MF) is used for vertical to slant TEC conversion and vice versa. There are several standards for the mapping function depending on the application (Figure 8). Here, we compare the MF recommended in MOPS 229D (see 419, 2008), GPS ICD (see GPS Joint Program Office 2004), and the modified single-layer model mapping (MSLM) given by Dach et al. [40] with ray-tracing results. The MF algorithms can be found in the given references.

Higher-Order Effects Estimation
The GPS L1 and L5 signals are traced through about 2000 electron density profiles (see left image of Figure 7), and higher-order ionospheric propagation terms including ray path bending effects are computed for single-frequency and dual-frequency ionosphere free linear combination. Figure 9 shows the variation of the second-and third-order residual terms, excess path length due to geometric bending, and excess TEC variation with the elevation angle.
Our investigation shows that the higher the peak electron density value, the higher the ray path bending effects. Additionally, the geometric bending effects (both excess path length and excess TEC) significantly increase with the decrease of the peak density height. It is worth pointing out that the geometric bending effects reduce to zero at 90 • elevation angle ( Figure 10).
For the second-order term estimation, we used a receiver location at -70 • N latitude and 120 • E longitude with a receiver-to-satellite azimuth angle set to 1 • . The second-order effect significantly varies with the user location and user-to-satellite azimuth angle at a certain elevation angle. For a GNSS user in the northern hemisphere, the magnitude of the second-order term is largest when the signal is received from a satellite in a southward direction. However, for a user in the southern hemisphere, the situation is reversed, e.g., the largest effect is observed when the signal is received from a satellite in a northward direction [12]. Therefore, the receiver location and azimuth angle are selected in such a way that we get the highest values for the second-order term estimation.
The large variations of ray-tracing results come from broader varieties of profile shapes. Apparently, the residuals due to ray path bending (both excess path length and excess TEC) are much higher than the second-and third-order residuals at low elevation angle (e.g., <15 • ). However, they decrease very rapidly with increasing elevation angle and disappear at zenith. Although the second-and third-order residuals are less than bending residuals at low elevation angles, they become larger than bending errors at higher elevation angles (>45 • ). The second-order residual does not reduce significantly with increasing elevation angle and therefore cannot be ignored even at zenith. Similarly, the third-order residual does not vanish at zenith. We found that the second-and third-order residuals are much less for ray-tracing computation compared to the values plotted in Figure 6 (blue dashed line), especially at low elevation angles.
shows smaller STEC values. With the increase of the elevation angle, each MF model shows better performance.

Higher-Order Effects Estimation
The GPS L1 and L5 signals are traced through about 2000 electron density profiles (see left image of Figure 7), and higher-order ionospheric propagation terms including ray path bending effects are computed for single-frequency and dual-frequency ionosphere free linear combination. Figure 9 shows the variation of the second-and third-order residual terms, excess path length due to geometric bending, and excess TEC variation with the elevation angle. Our investigation shows that the higher the peak electron density value, the higher the ray path bending effects. Additionally, the geometric bending effects (both excess path length and excess TEC) significantly increase with the decrease of the peak density height. It is worth pointing out that the geometric bending effects reduce to zero at 90° elevation angle ( Figure 10). For the second-order term estimation, we used a receiver location at -70° N latitude and 120° E longitude with a receiver-to-satellite azimuth angle set to 1°. The second-order effect significantly varies with the user location and user-to-satellite azimuth angle at a certain elevation angle. For a GNSS user in the northern hemisphere, the magnitude of the second-order term is largest when the signal is received from a satellite in a southward direction. However, for a user in the southern hemisphere, the situation is reversed, e.g., the largest effect is observed when the signal is received from a satellite in a northward direction [12]. Therefore, the receiver location and azimuth angle are selected in such a way that we get the highest values for the second-order term estimation. The large variations of ray-tracing results come from broader varieties of profile shapes. Apparently, the residuals due to ray path bending (both excess path length and excess TEC) are much higher than the second-and third-order residuals at low elevation angle (e.g., <15°). However, they decrease very rapidly with increasing elevation angle and disappear at zenith. Although the second-and third-order residuals are less than bending residuals at low elevation angles, they Simulation studies (see Figure 10a-d) show that the slant TEC varies from 360 to 1190 TECU. The corresponding second-order phase delay varies from 111 to 180 mm and 46 to 75 mm for the GPS L5 and L1 signals, respectively, whereas the residual error for the combined L1-L5 signal varies from 71 to 115 mm. The third-order phase delay varies from 6 to 43 mm and 2 to 13 mm for the GPS L5 and L1 signals, respectively. The third-order residual error for the combined signal varies from 10 to 72 mm. The excess path error varies from 0 to 496 mm and 0 to 154 mm for the GPS L5 and L1 signals, respectively whereas the excess path error for the combined signal varies from 0 to 278 mm. The excess TEC error for the combined GPS L1-L5 signal varies from 0 to 555 mm (see red plot of Figure 10a). Figure 6 shows the second-and third-order residuals obtained by approximation formulas. The reason for the high values is that we have assumed a radio wave propagation along the magnetic field line (e.g., Θ = 0), which gives excessive high values for the second-order residual for high TEC values at low elevation angles. The right panel of Figure 9 shows that the total dual-frequency residuals for the GPS L1-L5 combination can reach up to 1000 mm at very low elevation angles. The comparison between the developed overbounding model and the ray-tracing results shows that the 5.33σ model overbounds the raytracing results for elevation angles equal to or greater than 3 degrees.

Conclusions
We presented a model to overbound ionospheric residuals that remain after the application of the ionospheric-free linear combination. The model is designed to fit the needs of the MOPS for the next generation of dual-frequency SBAS using Galileo [19]. The final model (defined by Equation (12) and depicted in Figure 7) is elevation-dependent, easy implementable in the receiver software, and provides a conservative estimate of the worse-case residuals with a risk probability beyond 10 −7 . The design of the model is based on worst-case estimates of the geomagnetic induction and TEC. The latter has been derived through the analysis of a number of geomagnetic storm events and exploratory analysis of TEC data that have been observed during the events. Finally, since a number of approximations (e.g., the mapping function, estimates of the peak electron density) have been used, our model has been validated by ray-tracing results.
Finally, it should be mentioned that the presented model has been approved by members of the EUROCAE (European Organization for Civil Aviation Equipment) and is included in the new Galileo MOPS for Satellite-Based Augmentation System Airborne Equipment [18]. Future activities will aim to provide a less conservative model that may give more margins in the integrity estimation e.g., for more demanding operations than LPV 200 or CAT I autoland.