An advanced residual error model for tropospheric delay estimation

Global navigation satellite systems (GNSS) are widely used for safety-of-life positioning applications. Such applications require high integrity, availability, and continuity of the positioning service. Integrity is assessed by the definition of a protection level, which is an estimation of the maximum positioning error at extremely low probability levels. The emergence of multi-frequency civilian signals and the availability of satellite-based augmentation systems improve the modeling of ionospheric disturbances considerably. As a result, in many applications the tropospheric delay tends to become one of the limiting factors of positioning—especially at low elevation angles. The currently adopted integrity concepts employ a global constant to model the variance of the residual tropospheric delay error. We introduce a new approach to derive residual tropospheric delay error models using the extreme value analysis technique. Seventeen years of global numerical weather model fields are analyzed, and new residual error models are derived for some recently developed tropospheric delay models. Our approach provides models that consider both the geographical location and the seasonal variation of meteorological parameters. Our models are validated with a 17-year-long time series of zenith tropospheric delay estimates as provided by the International GNSS Service. The results show that the developed models are still conservative, while the maximal residual error of the tropospheric delay is still improved by 39–55%. This improvement yields higher service availability and continuity in safety-of-life applications of GNSS.


Introduction
In safety-of-life (SoL) navigation applications using Global Navigation Satellite Systems (GNSS), the main performance parameter of interest is integrity. A protection level is calculated that bounds the positioning error even at very small probability levels to assess the integrity (Zhu et al. 2018). The user predicts the protection level on the basis of error models of each of the contributing ranging error sources. To ensure that the protection level is not underestimated, each of these separate error sources needs to be modeled conservatively (Ober 2003). Our study focuses on conservatively modeling one of these error sources: the residual tropospheric errors.
GNSS utilizes ranging between the satellites and the receivers to determine the position of the user with respect to a geocentric reference system. Ranging is realized by the time-of-flight observation of the satellite signals. Redundant range observations in the positioning equations can be used to detect the non-nominal operation of satellites using the receiver autonomous integrity monitoring (RAIM) approach. However, the efficiency of the RAIM approach is limited in detecting frequently occurring faults (e.g., multipath in the urban environment) or biases affecting all of the satellites (El Mowafy et al. 2020). The satellite signals travel through the atmosphere and suffer considerable delays in the ionosphere and the troposphere. Lupsic and Takács (2019) showed that ionospheric effects can be corrected with an accuracy of approximately 1 meter, even for single-frequency users of satellite-based augmentation systems (SBAS). Tropospheric delays are usually taken into consideration using empirical models in absolute positioning, such as the ones detailed in Douša et al. (2018), Boehm et al. (2014), or Rózsa (2014). Jan (2010) showed that in case of processing dual-frequency pseudorange observations, the contribution of the tropospheric delays is higher than the ionospheric delays for satellites below 14°.
The performance of these models must be evaluated to ensure that safety-of-life users obtain reliable positions from their GNSS receivers. Collins and Langley (1998) published a model for the estimation of the maximum residual tropospheric delay error based on North American radiosonde observations and showed that the maximum residual error equals 60 cm. Van Graas and Zhu (2011) showed that even in differential positioning heavy rainfalls can cause a relative error of 30 cm over a 5 km baseline between a ground station and an aircraft.
The "standard" models used for residual tropospheric delay error given in the current quasi-standard (RTCA 2006) are considered very conservative (McGraw 2012). Although this is advantageous from a safety point of view, it has a negative effect on the availability and continuity of the positioning service.
The tropospheric delay model used today in SoL applications is recommended by the RTCA (founded as the Radio Technical Committee for Aeronautics) and is described in the company's Minimum Operational Performance Standards (RTCA MOPS). It has an associated maximum vertical error of 0.12 m in terms of standard deviation globally (RTCA 2006). Although (RTCA 2006) does not specify how this value was obtained, it agrees well with the value found by Collins and Langley (1998).
Van Leeuwen et al. (2004) studied the validity of this model in the Netherlands and also concluded that the model is highly conservative. Thus, there seems to be room for developing a less conservative model without compromising safety.
In the near future, more demanding applications are expected to arise, since multi-frequency and multi-constellation use of GNSS suffer from ionospheric delays less than today. Tropospheric effects cannot be eliminated by satellite signals using different frequencies. They need to be considered with empirical models. This creates a demand for more accurate tropospheric residual error modeling and ensures its importance in approximating integrity while maintaining sufficient system availability. Rózsa (2018) introduces a new method for the assessment of the tropospheric delay model performance using extreme value analysis. This method utilizes numerical weather model data and considers both the geographical and seasonal variations in the tropospheric behavior. Rózsa et al. (2017) provide a new overbounding model of the residual tropospheric delay for the RTCA model, resulting in a significant improvement in the calculated protection levels.
This study extends the results to include two newer tropospheric delay models: the ESA Galtropo (Krueger et al. 2004) and the Global Pressure and Temperature 2 wet (GPT2w, Boehm et al. 2014) models. The developed residual error models consider both the seasonal and the geographical variations of the estimation accuracy. The situation examined in this study always refers to single-point positioning with the tropospheric delay models used in blind mode. This implies that tropospheric delays are estimated without in situ meteorological data, using the model's built-in meteorological parameters.
In the next sections, we briefly review the relationship between GNSS integrity assessment and residual tropospheric error modeling. We summarize the applied methodology utilizing extreme value analysis on zenith tropospheric delay residuals, which are calculated using a ray tracing technique (Boehm and Schuh 2003) on global numerical weather models. We discuss the realization of the advanced residual tropospheric delay error models with three different complexity levels. Finally, the models are validated using an independent set of zenith tropospheric delays estimated at permanent sites of the International GNSS Service (IGS).

GNSS integrity and the role of tropospheric residual modeling
In GNSS, the integrity of the position is ensured using the concept of the protection levels. Protection levels provide an upper bound of the positioning error, which is only to be exceeded with a tiny probability (typically in the range of one out of every 10 7 -10 9 operations). When the computed protection level exceeds the acceptance level for the operation that the user wants to execute, the user cannot trust the position solution and will have to either abort the operation or resort to alternative sources of positioning information. Protection levels are calculated using models on the GNSS ranging error sources. Such a model must overbound the real residual error to ensure the correctness of the protection levels when the residual errors are translated into the positioning domain.
In RTCA (2006)  and long-term corrections, 2 i,UIRE is the model variance of the slant-range ionospheric delay residual error, 2 i,air is the variance of the airborne receiver error, and 2 i,tropo is the variance of the tropospheric delay residual.
The value of 2 i,tropo is to be calculated as: where TVE denotes the uncertainty of the vertical residual error of the tropospheric delay estimation with a constant value of 0.12 m and m i represents the mapping function that depends on the elevation angle i of satellite i . Note that TVE does not consider any climatic and seasonal variations over the globe, potentially leading to an overly conservative model in many geographical regions. In our approach, we omit the uncertainties regarding the mapping function coefficient. Although these can be substantial (Lagler et al. 2013), our study focuses on the integrity of the zenith tropospheric delay estimation only.
Combining (1) and (2), one can calculate the variance of the total residual error of range observations, which provides input for the estimation of the horizontal (HPL) and vertical protection levels (VPL) based on the actual geometry of the ranging observations (RTCA 2006).
In our study, we consider two recent tropospheric delay models, which consider the climatic variations of the fundamental meteorological parameters to increase performance with respect to that of the RTCA model. These models are briefly introduced in the next section.

Studied tropospheric delay models
Although the most common tropospheric delay models can operate in a so-called site augmented mode, where the user can supply in situ meteorological data, as mentioned before, this study focuses on the blind mode of operation, in which such data is not used and meteorological parameters are derived by the receiver as a function of time and location.
As described in Krueger et al. (2005), the ESA Galtropo blind model is based on the equations and constants of the previously developed TropGrid model (Krueger et al. 2004). However, it uses a different modeling approach. Its meteorological parameters were derived from the ERA15 statistical reanalysis data sets provided by the European Centre for Medium-Range Weather Forecasts (ECMWF).
The seasonal variation of each atmospheric quantity (surface water vapor pressure, water vapor lapse rate and mean temperature) is modeled. Furthermore, the diurnal variation of water vapor pressure is taken into consideration as well. The zenith hydrostatic delay (ZHD) is computed using the Saastamoinen model (Saastamoinen 1972, (2) 2 i,tropo = TVE ⋅ m i 1973), while the zenith wet delay (ZWD) is estimated using TropGrid's physical model, taking into account the modeled mean gravitational acceleration at the given latitude. The zenith total delay (ZTD) is converted to a slantpath delay using the Niell mapping function (Niell 1996). The GPT2w blind model provides the user with the pressure, the temperature and its lapse rate, the water vapor pressure and its decrease factor, and the weighted mean temperature as functions of the geographical location and the day of the year. Thus, both the annual and semiannual variations of these quantities are modeled over a global grid. The ZHD is computed using the Saastamoinen model, while the ZWD is estimated from the model described in Askne and Nordius (1987). The zenith delays are converted into slant delays by using the Vienna Mapping Function 1 (VMF1) given in Boehm and Schuh (2003).
In the remaining sections, we develop a less conservative model by considering the seasonal and geographical variations of the meteorological parameters, without losing the overbounding property. Such a model has the potential to increase the number of operations for which GNSS can be used without compromising safety.

Methodology
Although the applied methodology is discussed in detail in Rózsa (2018), it is briefly summarized here for the sake of clarity. The integrity of GNSS positioning services must be evaluated for extremely rare events with a significance level of 2·10 −7 according to the ICAO (2006) for various approach types. Assuming the duration of an average approach of 150 s without any concurrent approaches at the same time, the recurrence interval of an integrity event would be approximately 25 years.
Since no continuously available stationary error samples are available for the performance analysis using such a significance level, a probabilistic approach is used. Moreover, it must be considered that the residual tropospheric delay error usually does not follow the normal distribution at the tails of the distribution (Rózsa 2018).
To deal with this problem, we applied the generalized extreme value theory to characterize the tails of the distribution of the residual tropospheric delay error. This mathematical approach is widely used in the prediction of flood levels, but recently it has also been utilized for evaluating the integrity of GNSS and its augmentation systems (Ober et al. 2014;Rózsa et al. 2017;Rózsa 2018). In the next sections, more details are given on the mathematical formulation of our residual tropospheric delay error models.

Formulation of the residual tropospheric delay error model
To study the worst-case performance of tropospheric delay models, the empirical residual tropospheric delay model error values must be calculated first. The tropospheric delay estimates provided by the studied models are compared to "true" values. The resulting residuals are then analyzed using the generalized extreme value theory to formulate an overbounding model to estimate the maximal residual tropospheric error.
Since the "true" tropospheric delays are not known, numerical meteorological models are integrated to calculate these values in the zenith direction. Note that all our analysis is done for the zenith direction only and additional errors introduced by the mapping function are not taken into consideration.

Meteorological data
Various forms of meteorological data sets can be used for our purpose. Collins and Langley (1998) derived a residual error model for the UNB3 tropospheric delay model using North American radiosonde observations and showed that the maximum residual error is in the order of 60 cm for the zenith tropospheric delays. However, the geographical sparsity of the available radiosonde stations prevents us from deriving a reliable residual tropospheric error model that considers climatic variations for the whole planet. Therefore, global numerical weather models were chosen for our studies. While the superior resolution and accuracy of localized or regional weather models could yield more accurate reference values, these models only provide ZTD reference values with limited geographical validity. Hence, we did not use such local models. The most important requirements regarding the meteorological data used can be summarized as follows: • a sufficiently long time series of data sets must be available for the studies since the extreme value analysis incorporates an extrapolation of the distribution at the tails; • the meteorological data should be stationary and consistently processed in the study period; • the data sets must reproduce the seasonal and climatic variations of the meteorological parameters.
In order to fulfill these requirements, the ECMWF ERA-Interim global atmospheric reanalysis data sets (Dee et al. 2011) were used with a spatial resolution of 1° by 1° in latitude and longitude. Although model runs are available since 1979, our studies focused on a 17-year period between January 1, 2000 and December 31, 2016 using the four daily analyses referring to 00, 06, 12 and 18 UTC, respectively. Compared to Collins and Langley (1998), where the extrapolation to a recurrence time of 25 years was based on only 10 years of radiosonde observations, using data from the recent 17 years decreases the extrapolation effects of the extreme residual error and reduces the effects of climate change on our results.
To estimate the slant tropospheric delays in the study period, we ray-traced the numerical weather models between the receiver and the satellite. Thus, the ambient temperature, the dew point temperature and the geopotential values at all the available 37 pressure levels were obtained from the ECMWF. Afterward, the geopotential values were converted to altitudes, and the water vapor pressure values were calculated at each grid point.
The neutral atmosphere has a significant impact on the signal propagation even above the highest level of the ECMWF models. Rózsa (2014) showed that the hydrostatic delays of the atmosphere above 36 km altitude reaches the level of 1 cm. Thus, the ECMWF data were extended to the height of 86 kilometers using the International Standard Atmosphere parameters (ISO 1975). The heights and the corresponding meteorological data were also oversampled in a manner given in Rocken et al. (2001).
Finally, the obtained vertical atmospheric profiles of air pressure, water vapor pressure, and ambient temperature values were interpolated to a denser vertical profile to minimize the omission error of the numerical integration.

Residual error modeling
The methodology of the residual error modeling follows the method introduced in Rózsa (2018). To carry out the modeling, first, the real tropospheric delays were calculated using the ECMWF numerical weather model levels and the ray tracing technique (Boehm and Schuh 2003). The ray tracing was conducted for every grid point in the zenith direction for both the hydrostatic and the wet component of the tropospheric delay.
Then, the tropospheric delays obtained from the empirical models were subtracted from the reference values to obtain the residual error ( ). Due to seasonal effects, the residual errors of the tropospheric delays show a strong seasonal variation all over the globe. This is presented in Fig. 1, which depicts the time series of the hydrostatic delay residuals of the ESA Galtropo model for the latitude band between N41° and N50°.
To apply the extreme value theory for the estimation of maximum residual errors, we normalized the data set by an appropriate time-dependent function, which describes the seasonal variation of the standard deviation of the daily residuals. Since both the ESA Galtropo and the GPT2w models consider annual and semiannual variations of the meteorological parameters, a similar model was fitted to the variation of the daily standard deviations of the residuals (Fig. 2): In the equation above A 0 , A 1 , A 2 are the mean, the annual and the semiannual amplitudes of the daily standard deviations, DOY denotes the day-of-year, and DOY 0 represents the phase component. The normalization is carried out using a zero-mean assumption: In the equation above, corresponds to the residual error, while δ n denotes the normalized residual error. Although the zero-mean assumption is not entirely valid for the residuals, this simplification is necessary for formulating a residual error model that is consistent with the RTCA recommendations. Figure 2 shows the histogram and the normal probability plot of the normalized residuals. One can easily notice that the normalized hydrostatic delay residuals follow a heavy-tailed distribution especially for the negative residuals. The general extreme value (GEV) analysis technique provides a suitable method to extrapolate the extremes to any probability level for not only heavy-tailed distributions, but for light-tailed ones as well.
Thus, in the next step the GEV analysis is carried out on the normalized residuals. The annual extremes of the normalized residual error values were selected for both the minima and maxima for the GEV analysis, and GEV distributions were fitted to both the annual minima and annual maxima. Afterward, the positive and negative maximal residual errors representing the chosen recurrence time were estimated from the fitted distributions. Finally, the maximum expected error of the normalized residuals ( Δ n,max ) is defined as the highest absolute value of the positive and negative maximal residual error.
Protection levels are calculated from the standard deviation of the ranging errors using the law of error propagation and an appropriate safety coefficient (K) in GNSS integrity analysis calculations. Therefore, the maximum expected error is scaled to the standard deviation by: where σ n,max is the maximum standard deviation of the normalized residuals and K = 5.33 is the value of the inverse of the standard normal cumulative distribution function at the probability level of 1-10 −7 . Thus, the maximum standard deviation could be calculated in case of unbiased residual tropospheric delay error as the product of (3) and (5): It must not be forgotten that the normalization was done with the zero-mean assumption. In the case of biased residual tropospheric delay error, this assumption may lead failing to overbound the residual error. To avoid this, we increased the calculated standard deviation in (6) by a constant value depending on the maximal daily bias ( Δ 0 ) in the original residual error ( ). The motivation for adding this offset term is the following. Since a zero-mean assumption was used for the normalization of the residuals as well as the estimation of the seasonal variation of the maximal residual error, a biased residual error could exceed the maximum residual error estimated from the normalized residuals (δ n ). Therefore, (6) representing the seasonal variation of the maximum residual error must be increased by a bias term ( Δ 0 ∕K ). The value of Δ 0 is estimated by fitting another GEV distribution to the annual extremes of the daily mean values of the original residual error and estimating the maximum daily bias of the residual tropospheric delay error using this GEV distribution with the recurrence time of 25 years.
The maximum standard deviation is finally calculated as the sum of the maximum standard deviation computed from the normalized residuals (5), scaled by the factor (DOY) that captures the seasonal variations (3), and an offset ( Δ 0 ∕K ) that equals a scaled version of the maximum daily bias ( Δ 0 ) of the original residual error: The calculation is done as follows. First, the daily mean value of the residuals is calculated and the annual positive and negative extremes are selected over the study period. In the second step, a GEV distribution is fitted to the positive and negative residuals, respectively. Finally, the maximum bias of the residuals is estimated using the GEV distribution ( Δ 0 ) with 25 years of recurrence time to have a consistent result with the estimation of n,max and it is scaled to standard deviation using the same K coefficient as in (5). Finally, we obtain the standard deviation of the overbounding residual tropospheric delay error by combining the equations above: where the first part is the bias and the second part is the seasonal fluctuation multiplied by the scaled maximal expected value of the normalized residuals.
To develop residual tropospheric delay error models that can readily be incorporated in receiver firmware, we realized these models at three levels of complexity. In the next section, we study this and various other aspects of the implementation of the proposed approach.

Advanced residual tropospheric error models (ARTE)
Although the residual error modeling could be done for any specific grid point of the input numerical weather model data, it is advantageous to derive the simplest models possible, since they need to be implemented in receiver firmware. One of the advantages of the RTCA MOPS recommendation is its ultimate simplicity, as it provides just a single global value to model the zenith delay error.
To maintain the simplicity of the model to the extent possible while providing a less conservative but still safe model, three types of the ARTE model have been derived: in each 10° wide latitude band. It neglects the seasonal variation but does consider geographical variation and only requires the approximate latitude of the receiver as an input. • ARTE global constant model (ARTE-GCM): this model provides a global constant for the estimation of the maximal residual tropospheric delay error. Although this is the simplest model, naturally, it neglects both the seasonal and geographical variation of the residual error. If one uses the ARTE GCM model, no input data are needed to estimate the maximum residual error.
The three levels of ARTE models have been developed for a variety of tropospheric delay models. However, we focus on the results obtained for the ESA Galtropo and the GPT2w models. More information concerning the integrity analysis of the RTCA tropospheric delay model can be found in Rózsa et al. (2017).

ARTE band seasonal model
The ARTE-BSM model provides all six parameters ( Δ 0 , A 0 , A 1 , A 2 , DOY 0 , n,max ) in (8) for each latitude band of 10°. Tables 1, 2, 3, and 4 contain the parameters of each latitude band for the hydrostatic and wet components of the two studied models. The advantage of this model is that it considers both the geographical and the seasonal variations of the maximum residual tropospheric delay error, as depicted in Figs. 2 and 3. Moreover, it provides the residual error in the zenith direction for both the hydrostatic and the wet components separately. This opens up the possibility to use different hydrostatic and wet mapping functions for computing the residual slant delay for the respective component.
Looking at Figs. 3 and 4, we can see substantial differences between the seasonal variations of the different latitude bands. Moreover, in the case of the hydrostatic delays the level of the residual error is much lower for the tropical region than for the mid-latitude regions. Both phenomena can be explained by the annual variation of the air pressure in the different regions. In the tropics, the annual variation of the monthly mean air pressure is rather low and shows clearer annual and semiannual variations than the mid-latitude bands do. As an example, Fig. 5 depicts the surface pressure variation at a tropical site in Hong Kong and a continental site in Budapest. It can clearly be seen that the pressure variation pattern at the station in Budapest is more complex. Since hydrostatic delays can be calculated as a direct function of the observed or modeled air pressure, the hydrostatic delays show larger residual error in mid and high latitudes, especially in the wintertime, due to the more complex annual variation of the air pressure. This phenomenon can be observed on both hemispheres.
Concerning the wet component, the results show similar effects in the amplitude of the annual variations. The residual errors in the tropical regions show lower annual variations than mid and higher latitudes. However, the level of residual error is higher than it is for the higher latitudes during most of the year. This can be explained by the humidity Table 1 Parameters of the ARTE-BSM for the hydrostatic tropospheric delay computed using the ESA Galtropo model and its annual variations in the different regions. It can also be noted that the wet component shows smaller annual variations for the southern hemisphere. The reason is not clear. Possibly, this is caused by the lower accuracy of numerical weather models in the southern hemisphere (Simmons and Hollingsworth 2002) or even by the fact that the continental areas-subject to higher weather variations-are mainly located on the northern hemisphere.
The results also show that the hydrostatic and wet components of the tropospheric delays have similar magnitudes   Table 3 Parameters of the ARTE-BSM for the hydrostatic tropospheric delay computed using the GPT2w model of residual error. This might be surprising since it is known that hydrostatic delays can be modeled to millimeter-level accuracy when they are estimated using on-site air pressure observations. However, the empirical pressure functions used by the GPT2w and the ESA model are only capable of representing annual and semiannual variations. Although this representation is suitable for many regions of the world, it has deficiencies in the continental areas, since in these regions, air pressure variations are more complex. Moreover, residual error modeling estimates the magnitude of extreme errors, which is much higher than the standard deviation of the residuals, especially in case of heavy-tailed distributions (Fig. 2). This explains the fact that the estimated maximum residual error in the hydrostatic component has a magnitude similar to that of the wet component.

ARTE band constant model
Although the ARTE-BSM takes both geographical and seasonal variations of the uncertainties into consideration, the complexity of the model may cause some problems in its implementation in GNSS receivers. Therefore, simplified models were also developed. The ARTE Band Constant Model (ARTE-BCM) is derived from the ARTE-BSM by neglecting the seasonal variation and calculating the annual maximum of the uncertainties using (8) for each band i: Thus, 18 constant values for both the hydrostatic and the wet components must be stored in the receiver to calculate the maximal residual tropospheric delay error (Table 5). However, it must be noted that this model neglects the seasonal variation of the maximal residual tropospheric delay error. It is also important to realize that this simplification comes at the cost of higher protection levels. In certain latitude bands, this can potentially lead to reduced service availability in some seasons (Fig. 2, 3).

ARTE global constant model
The ARTE-GCM model is formulated by taking the maximum of the ARTE-BCM maximum residual tropospheric delay error values for both the hydrostatic and the wet components, respectively. According to Table 5, in terms of standard deviation, the maximal hydrostatic and wet residual tropospheric delay errors are 0.06 m and 0.07 m respectively. These values are lower than the RTCA recommendation (0.12 m) as the aggregated uncertainty is 0.092 m for the ARTE GCM model considering the residuals computed using the ESA Galtropo and the GPT2w delay models. If the same analysis is carried out for the RTCA tropospheric delay model, the final value agrees well with the RTCA recommendations (Rózsa et al. 2017).

Slant residual tropospheric error
According to the concept of the RTCA MOPS, the uncertainties must be calculated in the satellite direction. Thus, the vertical uncertainties calculated using the ARTE models must be converted to slant uncertainties. This can be done by using the appropriate mapping function recommended by the RTCA MOPS. Therefore, the maximal uncertainty of the slant total delay is: where m h and m w are the hydrostatic and wet mapping function values calculated for each respective satellite. Improving on the RTCA model, which only uses a single mapping function, the ARTE models could use separate mapping functions for the wet and dry components of the delay.
In the next section, we introduce the validation of the developed models using IGS zenith tropospheric delay data sets. The performance of the ARTE models is also examined. Both the conservatism of the model and the magnitude of the calculated tropospheric protection levels are evaluated. The former is a minimum requirement for safety-of-life applications, while the latter shows the improvement gained from better modeling the climatic effects.

Validation of the model performance using IGS reference data
In order to check the overbounding performance of the model, we validated the results with data acquired from the IGS (Dow et al. 2009). Our choice fell upon the data sets provided by the IGS as they can be considered independent from radiosonde and other sources that are incorporated into Hydr.

max (cm)
Wet max the numerical weather models. The data products containing total tropospheric delays are freely accessible through the IGS Global Data Center. Given that there are substantial differences in data quality between the stations, we chose 49 out of the initial 300, which have the least number of data gaps and are spread approximately evenly over the surface of the earth. The total tropospheric delay value at 12:00 UTC was used for the validation for each day from January 1, 2000 to October 31, 2017. Since there were no stations located between 81° and 90° of latitude north or south of the equator, the performance of these bands could not be validated.
The residual tropospheric delay error was computed by subtracting the estimated total tropospheric delay calculated using the ESA Galtropo and the GPT2w models from the values estimated at the IGS stations for each day in the aforementioned time period. These residual tropospheric delay error values were compared to the protection levels calculated using our three residual models ARTE-BSM, ARTE-BCM and ARTE-GCM for the respective tropospheric delay model.
To visualize the model performance, we created Stanford plots (Tossaint et al. 2006) for each latitude band and tropospheric delay model. As an example, we present two of these plots corresponding to the latitude band between N50°-N41° in Figs. 6 and 7. Approximately 25,000 data points of 5 IGS stations were used for the validation of the residual model performance in this given band. The figures show the relationship between the residual error computed using the IGS reference data (horizontal axis) and the protection level estimated by the ARTE-BSM (vertical axis) for the ESA Galtropo (Fig. 5) and the GPT2w (Fig. 6) tropospheric delay models, respectively. The plots demonstrate that the ARTE-BSM model provides better performance compared to the original RTCA recommendation, yielding no misleading data points, while lowering the current constant contribution to the protection level (0.64 m-tropospheric component only) by at least 39%.
To numerically describe the model performance in a concise manner, we calculated the ratio between the tropospheric residual error and the corresponding protection level for every IGS station involved in the validation. The maximum value of this ratio ( R max ) is presented in Table 6 for each latitude band. Table 6 shows that the developed residual models are successfully validated, and the highest ratios reach approximately 0.6, indicating that there is still a substantial amount of "reserve" contained in the calculated protection levels.
The results show a rather homogeneous performance over the globe, without any clear latitude dependency in the ratio between the residual error and the protection level. This suggests that the derived residual error models capture both the Fig. 6 Visualization of the ARTE-BSM model performance for the ESA Galtropo tropospheric delay model. For comparison, we also show the tropospheric protection level given by the RTCA. It is clear that there are no misleading data points where the computed protection level is lower is than the actual residual error seasonal and latitude-dependent variations of the residual errors. Compared to the current constant protection level defined in the RTCA recommendations, the ARTE-BSM model lowered the calculated protection levels by between 39 and 55% (depending on the latitude band) for both the ESA Galtropo and the GPT2w tropospheric delay models. This leads to higher availability and continuity in the positioning service.

Results and conclusions
We proposed a new tropospheric residual error model derived from ray tracing numerical weather model data over a period of 17 years. The resulting ARTE models were realized on three different complexity levels depending on whether the user chooses to take into account the seasonal and latitude dependence of the estimation accuracy. We have presented options to either neglect only the seasonal dependency or both the seasonal and the latitude dependency. We showed that in the latter case, the model is in good accordance with the de-facto integrity standard in aviation, the RTCA MOPS model. We also showed that the new model is much more efficient if seasonal and latitude dependencies are accounted for, which can lead to improved availability and continuity of the GNSS positioning service. Fig. 7 Visualization of the ARTE-BSM model performance for the GPT2w tropospheric delay model. For comparison, we also show the tropospheric protection level given by the RTCA Table 6 Maximal ratio between the hydrostatic tropospheric residual error in the zenith direction and the corresponding protection level given by the ARTE-BSM for both tropospheric delay models To evaluate the performance of the proposed models, we used residual errors computed using zenith total tropospheric delay estimates provided by the IGS in the validation process. These values were compared to the protection levels calculated by the ARTE models, and their overbounding characteristics were studied over a long period of time and throughout different latitude bands.
The proposed model proved to be more efficient than the currently used approach presented in the RTCA MOPS model providing smaller error estimates while maintaining sufficient overbounding of the residual errors as required for safety-of-life applications and other use cases demanding a high level of integrity.
It must also be noted that the proposed methodology can be applied to any other advanced tropospheric delay models as well.