WAAS and the Ionosphere – A Historical Perspective: Threat Model Evolution

2 SUPERTRUTH

Generating accurate ionospheric truth data is critical to the operation of satellite-based augmentation systems. Such truth data are required for the construction the undersampled ionospheric irregularity threat model and for continuous off-line analysis of ionospheric disturbances that might mandate its upgrade. The threat model has been updated six times over the lifespan of WAAS. The original threat model in IOC and the threat models in each of its updates have been derived from 1-s dual-frequency pseudorange and carrier-phase observations recorded at WRSs under disturbed ionospheric conditions. Stored in receiver independent exchange (RINEX) format (International GNSS Service, 2024), these measurements have been post-processed using the Jet Propulsion Laboratory (JPL) GNSS-Inferred Positioning System and Orbit Analysis Simulation Software (GIPSY-OASIS) package (Zumberge et al., 1997) and the Global Ionospheric Model (GIM) software package (Mannucci et al., 1998) to produce estimates of the total electron content (TEC) along GPS signal raypaths. The post-processing addresses four objectives (Komjathy et al., 2005): (1) data editing that includes identification of cycle slips, which subdivide each track into arcs, (2) leveling the carrier-phase measurements in each arc to the corresponding pseudorange measurements, (3) removal of satellite and receiver interfrequency hardware biases from relative TEC estimates, and (4) detection and elimination of spurious measurements using the observational redundancy provided by having three co-located receivers (or wide-area reference equipment) at each WAAS station. After the “best” estimate of the true TEC at each time step has been selected from the three threads of data, the resulting sets of estimated TEC data are designated as supertruth to distinguish them from truth data derived from measurements recorded by a single receiver at each WAAS station.

The identification of cycle slips on GPS L1-L2 phase observables and the removal of short arcs have been performed by a GIPSY data editor. Since supertruth is intended for ionospheric analysis under disturbed conditions, the GIPSY criterion used to identify each cycle slip has been set loosely to permit the retention of more data than would be kept using a criterion set for nominal conditions. Leveling accuracy depends strongly on arc length; thus, a tight cycle-slip criterion tends to degrade this accuracy, since the presence of many unwanted cycle slips leads to an excessive number of short TEC arcs. For similar reasons, arcs of less than 5-min duration have been dismissed. Cycle slips that go unflagged would, of course, corrupt the generated supertruth. Komjathy et al. (2005) note, however, that cycle slips are typically very large compared to the difference at adjacent points in a phase observable where no cycle slip is present; thus, this editing process should be efficient in identifying cycle slips. To mitigate multipath error on code measurements, a 5-s smoothing window is applied to the pseudorange observations.

The level of each continuous carrier-phase arc is established by averaging code minus phase ionospheric observables using elevation-dependent weighting. To determine the level of an arc prior to the availability of L2C and L5 signals, higher elevation data were first weighted more heavily using a model of P-code scatter that varied with elevation angle. If the actual P-code scatter computed from the measurements was found to deviate significantly from the leveling model, the corresponding portions of the data arc were discarded. (The noise characteristics of L5- and L2C-capable receivers could now be used to repeat the analysis, but this algorithmic change has not been deemed necessary.) A 5° elevation cutoff is applied to minimize the adverse impact of uncertainty in low-elevation measurements. For very short arcs, the elevation weighting has almost no effect on the weighted average.

After the observations have been edited, leveled, and decimated (i.e., downsampled from 1-s to 5-s intervals), satellite and receiver differential biases are removed from the carrier-phase observables to obtain unbiased line-of-sight TEC estimates. To determine these biases to high precision, the supertruth algorithm uses the thin-shell model of the ionosphere, constraining measurements on quiet days to satisfy the following observation equation (Mannucci et al., 1998):

STECrs=F(α,h)∑i=1NCiBi(θ,ϕ)+br+bs6

where F(α, h) is the thin-shell model obliquity factor that converts vertical TEC to slant TEC for raypath elevation angle α and shell height h (see Equation (4) of Sparks et al. (2022)), Bi(θ,ϕ) is one of N horizontal basis functions, each evaluated at the IPP solar geomagnetic latitude θ and longitude ϕ, C_i is the i-th basis function coefficient, and b_r and b_s are the receiver and satellite differential biases assumed to be constant over each 24-h period. (Here, geomagnetic coordinates are calculated according to a single iteration of the scheme of Bowring (1976).) GIM performs sequential least-squares parameter estimation by applying a Kalman filter to the carrier-phase-leveled ionospheric GPS observables and solving simultaneously for the basis function coefficients (C_i) and instrumental biases (b_r and b_s) for all three receivers at each WAAS station. The horizontal local basis functions B_i (θ, ϕ) employed in GIM are defined using a bicubic spline technique developed at JPL (Lawson, 2023). Each receiver at a WAAS station determines a distinct estimate of the TEC along the raypath connecting the station to the satellite. In this way, a redundant set of three TEC estimates is obtained for the TEC along any given raypath.

Modeling the distribution of charged particles in the ionosphere as a thin shell can serve as a significant source of error, even under the quietest and most nominal of ionospheric conditions. Two types of error afflict the accuracy of an estimate of the vertical TEC at an IPP when it is derived from a slant TEC measurement using the thin-shell obliquity factor (Sparks, 2013). First, the thin-shell model implies that the electron density is azimuthally symmetric near the IPP; distinct measurements that share a common IPP can generally produce compatible estimates of vertical TEC at the IPP only if the slant TEC associated with each raypath intersecting the IPP is independent of the raypath azimuthal angle at the IPP. Second, TEC error can result from a suboptimal choice of shell height. Since the thin-shell obliquity factor varies monotonically with shell height, the shell height that permits conversion of a given slant TEC measurement to the correct vertical TEC at the IPP is unique and is not generally known prior to the conversion. The magnitude of the error arising from these sources is small at mid-latitudes under nominal conditions. This tends not to be the case, however, at low latitudes (Rajagopal, 2004) or at mid-latitudes under disturbed conditions (Datta-Barua, 2004).

In the final post-processing step, the supertruth algorithm employs a voting scheme that compares the redundant observations from the three receivers on an epoch-by-epoch basis at 5-s intervals; the algorithm assigns a TEC value from one of the threads to the supertruth data set only if the measurements in the epoch from all three threads produce leveled ionospheric observables that agree to within a specified tolerance. When measurements in a given epoch from all three threads are present, the TEC values generated must agree to within 20% and 10 TEC units (TECU, where 1 TECU = 10¹⁶ electrons/m²); if measurements from only two threads are available, they must agree to within 40% and 20 TECU. TEC data that fall outside specified upper and lower TEC bounds are also ignored. Finally, a bound on the TEC scatter of leveled phase ionospheric observables is set loosely to 10 TECU to minimize data loss. The minimal error that remains at the conclusion of the computation is due primarily to phase multipath, receiver noise, uncorrected cycle slips, and uncertainties in the estimation of receiver and satellite biases.

To facilitate data access, each ASCII supertruth data file is converted into a binary MATLAB file prior to constructing a threat model. An additional decimation is generally performed at the time of this conversion. The user also has the option of creating a threat model using a subset of the data in the MATLAB file by specifying a “step” value indicating the number of seconds between epochs at which fits are performed.

To date, five sets of WAAS supertruth files have been generated, and each set has been assigned a distinct version number. Table A1 specifies the version number for the observational data file set used to construct each released threat model. Version 1 supertruth consists of all the TEC data files used to construct the IOC ionospheric threat model, spanning the period January 11, 2000, to March 31, 2001.

Version 2 supertruth extended the period of coverage to include storm data through the end of solar cycle 23; the final storm day in this set is January 22, 2005. (In Version 2 supertruth, all Version 1 data files were recomputed except that of January 11, 2000.) Version 2 supertruth incorporated various algorithmic improvements, achieving, on average, a 30% increase in the volume of data retained in a 24-h period over that of the previous algorithm (Komjathy et al., 2005). Of perhaps greater significance, the Version 2 algorithm allowed critical data to be recovered from periods of significant ionospheric disturbance. Furthermore, the Version 2 algorithm improved TEC accuracy by extending the period of the data processed by 6 h both prior to the start of the day and after the day’s conclusion. This permitted an improvement in the leveling of the data at the day boundaries. The threat model of LPV Release 6/7 was based on Version 2 supertruth.

Version 3 supertruth consists of TEC data from only six dates: July 26–27, 2004; November 8–10, 2004; and January 22, 2005. The need to recompute the supertruth associated with these dates arose after the discovery of a minor software error that had caused data to be omitted from the corresponding Version 2 supertruth files on these dates after a new satellite was added to the GPS constellation. The recomputation caused slight changes in the retrieved hardware biases, which in turn, caused a small modification of the values of the fit residuals used to construct the ionospheric threat model. Version 3 supertruth was first implemented in the generation of the threat model assigned to LPV Release 8/9.2 in September 2008. The set of supertruth files used to generate this threat model was retained in the tabulation of the threat models for WAAS Follow-On (WFO) Release 3 and Release 46-CY16.

A major modification of the supertruth algorithm occurred in 2011 when data-editing parameters were altered so that additional data, previously eliminated by the editing performed by JPL’s GIPSY program, were included in the supertruth file for each storm day. When the updated algorithm was applied to observational data from solar cycle 23, the resulting output was designated Version 4 supertruth. With the advent of solar cycle 24, yet another algorithmic modification was implemented: Version 5 supertruth employed substantially the same algorithm as Version 4 supertruth except that the period of the data processed both prior to the start of the day and after the day’s conclusion was extended from 6 h to 9 h. An ionospheric threat model constructed using Version 4 and Version 5 supertruth was not fielded until Release 51-CY18.

3 MEXICAN DATA

The computation of the IOC ionospheric threat model was based entirely upon historical observations recorded by receivers at the 25 WRSs that comprised the WAAS network prior to WAAS’s commissioning in 2003. In 2007, LPV Release 6/7 expanded the network at lower geomagnetic latitude to include stations in Mexico. Since the ionosphere exhibits greater variability at low latitude than it does at mid-latitude, it became necessary to acquire historical data from Mexican stations to ensure that a threat model constructed from data recorded at mid-latitude continued to protect users at low latitude. With assistance from Patricia Doherty of Boston College and from various Mexican institutions, notably the Universidad Nacional Autόnoma de México and the Instituto Nacional de Estadística y Geografía, ionospheric delay measurements observed on storm days of interest were collected from various dual-frequency Mexican receivers and integrated into the corresponding supertruth data.

Unlike each WAAS station, only one receiver resided at each Mexican station. Therefore, these Mexican data were considered truth data as opposed to supertruth data derived from observations of three co-located receivers. The solar-cycle-23 WAAS supertruth integrated with Mexican truth was designated a WAAS+Mexico data set to distinguish it from the WAAS-only data set consisting of the same data without the Mexican contribution. A consequence of variability in the quality of the Mexican measurements was that data from a different subset of Mexican receivers were adopted for each storm day. Historical sites near current WAAS WRSs were chosen preferentially to mimic the observations that a Mexican WAAS receiver might have recorded had it been present on a given storm day. (The quality of the Mexican data and the actual sites selected for each storm day have been discussed by Paredes et al. (2008)). The Mexican data were preserved at epochs separated by 30 s; in the integrated supertruth data files, the WAAS-only data were decimated to match the decimation in the Mexican data.

The presence of receivers in Mexico complicated the process of obtaining accurate vertical delay estimates at southern IGPs. Especially when viewing the southern sky, these receivers recorded measurements whose IPPs often lay much further south than the southernmost IPPs included in IOC fits. IPPs occurring at low geomagnetic latitude occupy a geographic region where the ionosphere is typically much more structured than at mid-latitude, and the assumed planar model of ionospheric delay is no longer valid. To maintain estimation accuracy, it is necessary to eliminate such measurements from the fits at local IGPs. In LPV Release 6/7, this elimination was accomplished in two steps. First, the receivers in question were identified and placed in a low-latitude WRS mask. Second, the nominal border of the region of high ionospheric variability was specified by two connected line segments, initially designated the Mexican deprivation line and now called the pierce point filter line (PPFL), defined by endpoints whose eastern, middle, and western latitudes and longitudes are labeled λ_E, ϕ_E, λ_M, ϕ_M, λ_W, and ϕ_W, respectively. The first of these two segments closely followed local contours of geomagnetic latitude. For measurements observed by the masked receivers, those whose IPPs lay south of the PPFL were excluded from all fits. This methodology has continued to be used when including solar-cycle-23 data in the threat models released after that of LPV Release 6/7. In subsequent threat models, however, the position of the PPFL has been adjusted. In the first companion paper (Sparks et al., 2022), the position of the PPFL for each threat model release is plotted in Figure 3, and coordinates describing the position of the PPFL are listed in Table A1.

4 USER SMOOTHING

Due to post-processing, supertruth data may be regarded as providing a more accurate estimate of the true slant delay at the GPS L1 frequency (1575.42 MHz) than the observational data determined directly by a WAAS user in real time. The latter is a delay computed using a phase-corrected lag filter. The equations describing this filter are as follows:

Pproj(ti)=Psmooth(ti−1)+λ2π(ϕ(ti)−ϕ(ti−1))7

Psmooth(ti)=αUSρ(ti)+(1−αUS)Pproj(ti)8

where ϕ(t_i) is a carrier-phase measurement (in radians) recorded at time t_i, ρ(t_i) is the corresponding pseudorange measurement (in meters) recorded at time t_i, λ is the wavelength (in meters), and α_US is a smoothing constant. P_smooth (t_i) is the smoothed pseudorange (in meters) adopted by the user in a position computation at time t_i. The Minimum Operational Performance Standards (MOPS) for Global Positioning System/Satellite-Based Augmentation System Airborne Equipment (RTCA, 2016) specifies α_US as ∆d/100 to invoke 100-s smoothing with a data decimation interval of ∆d. Such a filter reacts more slowly to rapid ionospheric fluctuations that are captured in supertruth. In this sense, supertruth does not accurately reflect the user’s experience. For this reason, pseudorange-smoothing by the user’s equipment potentially represents an additional source of positioning error, a concern that has been designated the user smoothing threat.

LPV Release 6/7 was the first WAAS release to take the user smoothing threat into account when constructing the ionospheric threat model. To simulate the impact of such a threat, an estimate of the user-smoothed slant delay STEC_smooth evaluated according to Equations (7) and (8) at time t_i is first generated as follows:

STECsmooth(ti)=αUSSTEC(ti)+(1−αUS)(STECsmooth(ti−1)−STEC(ti)+STEC(ti−1))9

where STEC(t_i) is the slant delay (in meters at L1) along a track at time t_i as read from a supertruth file. In IOC, each residual needed to evaluate σ¯undersampled,κ2, as defined by Equation (3), was calculated twice: once using unsmoothed supertruth data and a second time using Equation (9). The larger of the two residuals was then incorporated into the tabulation of the threat model. This approach to treating the user smoothing threat has been retained in all subsequent releases; the value of α_US has not been modified.

5 SIMULATION OF IONOSPHERIC IRREGULARITY UNDERSAMPLING

Ionospheric irregularities can induce dramatic changes in TEC magnitude over comparatively short distances. Changes of greater than 20 m of vertical delay at the GPS L1 frequency have been observed by WAAS over a few hundred kilometers (Walter et al., 2004). A violation of the inequality in Equation (2) may occur when a user measures a GPS signal passing through a disturbed region of the ionosphere that is not sampled (or is poorly sampled) by signals recorded at WAAS receivers. To construct each branch of the ionospheric threat model, it is necessary to simulate such situations. This simulation is performed by treating each WAAS observation in a specified data set as a user measurement, estimating the corresponding TEC value using a fit centered on the IGP closest to the measurement IPP, and then evaluating the magnitude of the vertical delay difference specified in Equation (2). Fit residual errors are tabulated in this fashion for each branch of the threat model, with the goal of determining the maximum values that σ¯undersampled,κ2 can assume when the various detectors specified by the branch have not tripped.

Since distributions of fit IPPs associated with actual observations cannot be assumed to encompass the worst possible configurations of ionospheric threats, a technique designated data deprivation has been employed to simulate many additional data sets. This approach allows WAAS to capture better the effects of undersampling. To evaluate σ¯undersampled,κ2, fit residuals have been calculated not only for sets of vertical delays drawn from the entire set of slant delays actually measured by WAAS in each epoch but also for sets of vertical delays that systematically exclude specified measurements in the epoch so as to generate alternative samplings of actual ionospheric threats affecting the measurements.

Data deprivation has focused on two distinct types of undersampled threats: irregularities that reside in the interior of the coverage region and those located at the edge of coverage. Bearing in mind that the locations of WAAS WRSs are typically separated by hundreds of kilometers, the first type of threat can escape detection by being sufficiently compact; the second type may go unobserved by being situated largely outside the sampling region. Studies have shown that the latter type presents the more formidable challenge to WAAS integrity (Blanch et al., 2002; Sparks & Altshuler, 2021), and indeed, this type of threat generally provides the dominant contributions to the threat model.

For each IPP present in a specified epoch, the algorithm implementing data deprivation proceeds as follows: a specified set of IPPs are excluded from consideration, and fit IPPs in a domain centered on the IGP nearest the test IPP (i.e., the IGP in whose threat domain the test IPP is located) are selected from the IPPs that remain. The process of fitting the vertical delay at these fit IPPs to a plane is then repeated in the same manner that a fit estimates the IGD to be broadcast at the IGP. When the resulting σ¯undersampled,κ2 at the test IPP is positive, it is included in the tabulation of the threat model provided that the irregularity detector at the IGP has not been triggered, i.e., provided that the value of the irregularity metric based on the chi-squared statistic of the fit does not exceed a specified irregularity threshold.

In IOC, fit IPPs were removed simply according to their geographic location. Two data deprivation schemes, designated annular and three-quadrant, addressed, respectively, the interior and edge-of-coverage threats. In the first case, measurements were excluded from fits when their IPPs lay in a 200-km-thick annular mask centered on the IGP (see Figure 1(a)). For this set of deprivation masks, the data in each fit epoch were processed 10 times, beginning with a circular data mask extending from 0 to 200 km and then moving the annular domain successively outward 200 km at a time until the outer domain boundary reached 2000 km. This scheme was designed to account for the impact of irregularities spatially localized within the fit domain.

In the second scheme, measurements were excluded from fits when their IPPs were located in any of three quadrants of a rectangular (latitude–longitude) grid whose origin lay on a diagonal passing through the fit IGP in question. Figure 1(b) displays the mask boundaries defined when fits were performed using the measurement IPPs restricted to the southwest quadrant. The origin of the grid was positioned at each of 11 uniformly spaced locations along the southwest-northeast diagonal within a 10° × 10° rectangular domain centered on the IGP. Similar sets of masks were defined for fits performed using IPPs located only in the northwest quadrant, the northeast quadrant, and the southeast quadrant. As a result, the data of each fit epoch were processed according to this scheme 44 times. The three-quadrant scheme attempted to account for the influence of spatially extended irregularities that remain undetected due to a poor distribution of IPPs. This scheme was devised to reflect ionospheric behavior at the edge of the WAAS grid, particularly when disturbances approach the fit domain from outside the region of coverage.

The data deprivation masks implemented in IOC were later deemed to suffer from a serious limitation: the resulting IPP distributions used in the data-deprived fits were not representative of actual IPP distributions. For example, consider the three-quadrant mask depicted in Figure 1(b). Removing all IPPs in the yellow region produced a very sharp cutoff in the geographic extent of the IPPs included in the fit. At the edge of WAAS coverage, however, the actual distribution of IPPs present in each epoch tapers off more gradually. Thus, a large gradient moving from the yellow region into the southwest quadrant would be sampled in a manner that is not modeled accurately by a three-quadrant mask.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 1

Data deprivation masks used in IOC WAAS: (a) annular masks, (b) three-quadrant masks that exclude from the fit the IPPs located outside the southwest quadrant The red dashed lines define the geographic boundaries of different masks. The yellow region in each plot identifies where IPPs are excluded from the fit for one specified mask. Red dots indicate the locations of GPS signal IPPs. Green dots at the corners of grid cells represent IGPs.

In September 2007, LPV Release 6/7 introduced a new set of data deprivation schemes that sought to remedy this problem. Rather than masking fit IPPs according to their geographic location, IPPs were removed according to the WAAS stations that recorded their measurement. To address undersampled threats localized within the fit domain, LPV Release 6/7 replaced the annular deprivation masks of IOC with single-station deprivation – masks that each eliminated IPPs associated with measurements recorded by a single specified WAAS WRS. Edge-of-coverage threats were examined using directional station deprivation, i.e., by systematically removing fit IPPs generated by measurements recorded at the N_station most extreme locations in a specified geographic direction within the WAAS coverage region. In LPV Release 6/7, N_station was varied between 2 and half the total number of stations in the network. A sequence of such masks for removing IPPs was defined in each of the eight cardinal directions: north, south, east, west, northeast, southeast, northwest, and southwest. To determine the order of the stations to be removed, station locations were sorted according to (1) latitude for a north–south orientation, (2) longitude for an east–west orientation; (3) the sum of latitude and longitude for a northeast–southwest orientation, and (4) the difference between latitude and longitude for a northwest-southeast orientation. (In all fits, data recorded by the station at Honolulu, Hawaii, were ignored owing to the geographic isolation of Honolulu from all other WAAS stations.)

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Directional station deprivation in a network configuration consisting of the 25 IOC WAAS stations (i.e., the first 25 stations listed in Table A1 of Sparks et al. (2026))

The depicted mask removes from fits of vertical delay the measurements recorded by receivers at the 12 most northeastern WAAS stations (in red). Data recorded at the station in Honolulu (black triangle) are ignored as well.

Directional station deprivation allows us to simulate the impact of a threat located at the edge of coverage using measurements of signals penetrating a disturbance that in fact occurred in a well-sampled region of the ionosphere. Figure 2 displays an example of a station mask that implements directional station deprivation in the IOC configuration of WAAS stations. This mask removes from fits of vertical delay at all IGP the measurements recorded by the 12 most northeastern network stations (depicted in red). This station mask has given rise to one of the most influential critical points in the disturbed-time branch of the ionospheric threat model, a critical point generated from measurements recorded during the extreme storm of November 20–21, 2003, that causes a sharp jump in the magnitude of the tabulated overbound at R_fit = 1450 and RCM = 0.575 (the critical points of a threat model branch are defined in Section 6 of Sparks et al. (2026)).

Figure 3(a) shows the distribution of IPPs measured in an epoch during which WAAS detected a localized ionospheric irregularity just west of the Great Lakes region. An X marks the location of an IGP whose threat domain contained IPPs (pictured in Figure 3(a) but not Figure 3(b)) associated with signals that sampled this irregularity. Note that the distribution of fit IPPs tapers off gradually in Figure 3(b) as one moves northeasterly into the masked region. If a three-quadrant mask were used to mask an analogous region, the transition to the area of masked IPPs would be unrealistically abrupt. Removing fit IPPs according to station rather than geographic location improves the verisimilitude of the simulated IPP distributions.

In the absence of data deprivation, the irregularity visible in Figure 3(a) does not contribute to the ionospheric threat model. When applied to this epoch’s set of IPPs, however, the station mask displayed in Figure 2 produces a large value of σ¯undersampled,κ2 for at least one of the IPPs in the threat domain of the IGP in question. When the disturbed-time branch of the WFO Release 3A ionospheric threat model is tabulated, this σ¯undersampled,κ2 value ultimately proves to be the largest value encountered in the bin whose lower-left-hand corner lies at R_fit = 1450 and RCM = 0.575 (see Figure 12(b)).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Distribution of fit IPPs in the epoch on November 20, 2003, that generates a critical point at R_fit = 1450 and RCM = 0.575 in the tabulation of the disturbed-time branch of the ionospheric threat model: (a) all measurement IPPs in the epoch, (b) the IPPs that remain in the domain after directional station deprivation has been applied to remove the 12 most northeastern stations

In each plot, X marks the location of the IGP in whose threat domain (indicated in red) the IPP responsible for the critical point (pictured in Figure 3(a) but not Figure 3(b)) is located. The fit domain centered on the IGP is indicated by the dashed line in Figure 3(b).

It is clear from Figure 3(b), why directional data deprivation enables such a large value of σ¯undersampled,κ2 to form: none of the measurements associated with the fit IPPs sample the threatening irregularity, a storm-enhanced density plume (Foster et al., 2005). An irregularity, similar to the one represented in Figure 3(a), that lay a comparable distance off the actual coast of the coverage region would represent a serious threat to aircraft navigation safety. This is an irregularity the threat model must take into account. Note: the data deprivation mask used to generate Figure 3(b) prevents the ESD from having tripped before the epoch in question (Sparks et al., 2022). When all IPPs are present, however, the ESD trips prior to this epoch, removing the possibility of any threats from this epoch contributing to the ionospheric threat model.

The LPV Release 6/7 implementation of both single-station deprivation and directional station deprivation has been retained in the construction of all subsequent threat models. The maximum N_station in the latter deprivation scheme, however, has grown as the number of network stations has been increased: at present, the maximum N_station has risen from 12 to 19.

Each of the data deprivation schemes discussed above can be considered optimistic in the sense that the probability is small an undersampled threat of interest will align itself in the worst possible orientation with respect to a hole in the IPP coverage. Walter et al. (2004) proposed a more aggressive approach to data deprivation, a scheme tailored to irregularities that are located in the interior of the coverage region. The idea was to target for removal IPPs corresponding to signals that sample the ionospheric threat in question. The objective of this scheme was to remove only those points that give rise to the largest disparity from a planar fit, points that would otherwise alert WAAS to the threat’s presence. In this sense, the points removed may be regarded as outliers. As such points are removed, the remaining IPPs become more consistent with a retrieved planar fit characteristic of quiet conditions, and the likelihood that the threat will remain undetected increases. Since the relationship between the undersampled threat and the fit IPPs removed is no longer random, this scheme is designated malicious deprivation. In LPV Release 6/7, malicious deprivation was implemented in the data deprivation schemes that simulate holes in coverage, i.e., single-station deprivation. (Malicious deprivation was not applied to directional station deprivation; since directional station deprivation addresses threats posed by ionospheric irregularities at the edge of coverage, and malicious deprivation is relevant only to irregularities occurring at holes in coverage.) Malicious deprivation has been applied throughout the coverage region for deprivation masks defined on a station-by-station basis. Both one and two IPPs have been removed maliciously for each mask associated with a single station (two-point deprivation has not required that the points be physically contiguous).

It should be noted that, since the fit IPP selection algorithm is invoked again after a fit IPP has been identified for malicious removal, the revised fit at the IGP is generally performed with a new fit IPP added to restore the targeted number of IPPs. For this reason, the χirreg2 irregularity metric does not necessarily decrease monotonically with respect to malicious deprivation. As a conservative measure, therefore, the smallest χirreg2 metric from among the prior fit(s) and the current fit at each IGP in each epoch is chosen as the metric for determining whether the revised estimate of vertical delay at the IPP should be tabulated in the threat model. Since malicious deprivation used 24 stations in LPV Release 6/7 and two points in each fit were successively selected for malicious deprivation, 72 single-station deprivation masks were used to prepare LPV Release 6/7.

The implementation of single-station and directional station data deprivation in LPV Release 6/7 coincided with the introduction of the ESD into WAAS. The logic governing the triggering of the ESD was applied to each data deprivation mask independently: all IPPs associated with the WAAS stations excluded by a given mask were removed from the vertical delay estimation (and from the evaluation of the irregularity metric) at each given IGP for the entire run. Directional station deprivation effectively enabled a storm gradient that would otherwise trigger the ESD to be moved to the edge of coverage and remain hidden, thereby generating a “worst-case” threat. As discussed by Sparks et al. (2022), the ESD trip threshold T_ESD,trip, the ESD onset confirm interval t_ESD,confirm, the ESD recovery threshold T_ESD,recovery, and the ESD recovery confirm interval t_ESD,recovery were initially set to 50, 1 h, 42, and 8 h, respectively.

When an ionospheric disturbance causes the irregularity detector at a given IGP to trip, WAAS sets the GIVE broadcast at that IGP to GIVE_MAX. Random statistical fluctuations can cause the irregularity metric at that IGP to drop below the trip threshold while the local ionosphere remains in a disturbed state. To prevent a premature decrease in GIVE, WAAS continues to broadcast a value of GIVE_MAX at the IGP for a period designated the hysteresis time t_irreg,hyst. In LPV Release 6/7, this hysteresis interval was extended from its value of 15 min in IOC to 30 min. As a conservative measure, however, the hysteresis interval t_irreg,threat used to tabulate values of σ¯undersampled,κ2 in the ionospheric threat model was set in LPV Release 6/7 and in all subsequent releases to 15 min, i.e., 0.5 times t_irreg,hyst. (In IOC WAAS, t_irreg,threat was zero, i.e., there was no hysteresis period over which a triggered irregularity detector excluded threats from the threat model.) Thus, the irregularity detector is less aggressive at removing threats in the construction of the threat model than it is in the system. An additional conservative step from LPV Release 6/7 to the present has been to set the measurement noise covariance to zero in the construction of the threat model. (The methodology used to generate supertruth is designed to make measurement noise negligible.) Consequently, R_noise is reduced to 1, rendering the ESD and each irregularity detector in the threat model less likely to trip than in the system (Sparks et al., 2022).

The implementation of data deprivation in LPV Release 6/7 (including specification of a zero-measurement noise covariance and the given values of t_irreg,hyst and t_irreg,threat) has been preserved in all subsequent WAAS ionospheric threat models. WAAS assumes that directional station deprivation and single-station deprivation augmented by malicious deprivation, when applied to a historical storm data set used to generate an ionospheric threat model, are sufficient to model the worst behavior of the ionosphere that does not trigger a local irregularity detector or the ESD. When the MSD was later introduced, the same logic for handling malicious deprivation in the presence of the ESD was applied in the presence of the MSD.

6 FIT IPP DISTRIBUTION METRICS PARAMETERIZING THE IONOSPHERIC THREAT MODEL

As discussed by Sparks et al. (2022), each fit of vertical delay is conducted within a fit domain consisting of a region on the ionospheric shell enclosed by a circle whose points lie a Euclidean (straight-line) distance R_fit from the IGP (see Figure 4). To select measurements to be included in a fit, the WAAS search algorithm seeks at least N_target IPPs within a distance R_min of the IGP at the fit center. If fewer IPPs lie within this fit domain, the fit radius R_fit is increased until the circle it defines encompasses N_target points. If this circle fails to enclose N_target points when R_fit attains a maximum value of R_max, the fit is conducted using fewer than N_target points as long as the number of fit points exceeds a specified minimum N_min (otherwise, the GIVE monitor sets the broadcast GIVE at the IGP to “not monitored”). The process by which the values initially assigned to these parameters (R_min = 800 km, R_max = 2100 km, N_target = 30, N_min = 10) were selected is described by Sparks et al. (2026). These values have been retained throughout the subsequent evolution of the ionospheric threat model.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Schematic diagram illustrating how vertical delay is estimated throughout the construction of the ionospheric threat model

A delay estimate is evaluated at each IPP that lies within the threat domain surrounding a given IGP. Estimates are achieved by fitting all of the measured values of slant delay, converted to vertical delay, whose IPPs are located in a fit domain defined by a circle of radius R_fit centered on the IGP.

R_fit may be regarded as a proxy for the mean IPP density, which necessarily decreases as the fit domain expands to encircle the targeted number N_target of fit IPPs. For this reason, it was the first parameter proposed as a fit IPP distribution metric: as R_fit increases, progressively larger irregularities can go undetected.

The introduction noted that a threat to WAAS integrity is considered to exist whenever the ESD and the local irregularity detector have not tripped and the inequality in Equation (2) is violated. Prior to the commissioning of WAAS in 2003, the WAAS Integrity Performance Panel considered making the augmentation of the GIVE at an IGP simply a function of the magnitude of R_fit: a larger fit radius corresponds to a larger adjustment provided by the undersampled ionospheric irregularity threat model. For the storm data set analyzed in the construction of the IOC threat model, Figure 5 shows the number of ionospheric threats as a function of the fit radius R_fit and the magnitude of the vertical delay fit residual associated with each threat. Figure 5(a) shows the results of tabulating only those threats that arise due to actual measurements. In this case, the vast majority of actual threats occur when R_fit is near the specified minimum R_min = 800 km. Figure 5(b) shows the results obtained by adding to the tabulation of Figure 5(a) the results of simulating threats by applying the IOC data deprivation schemes, i.e., annular and three-quadrant masks, to the same set of observations. Note that on both the vertical scale and the color bar scale, the maximum values are much larger than in Figure 5(a). The final column includes each simulated threat for which the corresponding fit at the IGP incorporates fewer than N_target fit points. The spike in this column reflects the fact that the probability of large irregularities going undetected increases dramatically as the number of points included in the fit drops below N_target. An advantage of choosing R_fit as an IPP distribution metric is that when fit residual errors are binned as a function of R_fit, large residuals derived from fits of sparse measurements incorporating fewer than N_target fit IPPs are relegated to the final bin column, which is rarely accessed in the interior of the WAAS service volume.

Note, however, that according to Figure 5(b), the maximum residual delay error can exceed 4 m (if the fits of simulated IPP distributions are taken as representative of possible actual threats) even when the fit radius is as small as 1000 km (where at least N_target IPPs are always included in the fit). To characterize the worst-case user error by adopting an overbound of these tabulated data proves too conservative for successful WAAS operation. Basing the augmentation of the GIVE on the worst-case errors parameterized by R_fit alone results in WAAS always broadcasting large error bounds and, consequently, failing to achieve satisfactory availability.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

Ionospheric threats generated during storms included in the IOC threat model (i.e., delay fit residuals |I¯−I∼| that violate the inequality in Equation (2)), binned as a function of fit radius R_fit and the delay fit residual: (a) threat counts without data deprivation, (b) threat counts with data deprivation

To remedy this situation, IOC WAAS implemented a second IPP distribution metric designated the relative centroid metric (RCM), as defined in Section 1. When fit IPPs are distributed uniformly throughout the fit domain, this metric is nearly zero; the metric approaches unity when the fit IPPs congregate near a single point at the edge of the fit domain. Thus, the RCM may be taken as a measure of the degree of uniformity in the distribution of IPPs throughout the fit domain: a larger RCM value corresponds to a greater likelihood that an irregularity will be undersampled.

Tabulating ionospheric threats according to the relative centroid metric allows us to distinguish threats that arise when the distribution of fit IPPs within the fit domain is optimal from those that occur when this distribution is highly non-uniform and therefore prone to undersampling irregularities of larger magnitude. Figure 6 shows how binning data as functions of both the fit radius and the relative centroid metric enables us to differentiate threats within each bin plotted in Figure 5: Figures 6(a) and 6(b) present the two-metric analog of Figure 5(a), and Figures 6(c) and 6(d) present the same for Figure 5(b). Note that, the color bar scale in Figure 6(a) differs from that in Figure 6(c), and the color bar scale in Figure 6(d) saturates at 10 m, causing all pixels exceeding 10 m to be colored brown. Furthermore, note that the benefit of introducing a second IPP distribution metric (at least in IOC WAAS) is primarily a consequence of the impact of data deprivation on threat assessment. There would be little or no benefit to using the RCM if the magnitudes of the fit residuals binned for the simulated threats (included in Figure 6(d)) were distributed as evenly as they are for the actual threats alone (Figure 6(b)). Figure 6(a) shows, however, that under nominal conditions, the RCM is typically found to be less than 0.4 (except possibly at IGPs near the edge of the WAAS coverage). Thus, using the RCM as a fit IPP distribution metric reduces the GIVE magnitude that would otherwise be required to protect the user from undetected threats occurring even when the RCM is less than 0.4.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

Tabulation of ionospheric threats generated during storms included in the IOC threat model (i.e., delay fit residuals |I¯−I∼| that violate the inequality in Equation (2)), binned as a function of fit radius R_fit and RCM: (a) threat counts without data deprivation, (b) maximum delay fit residuals without data deprivation, (c) threat counts with data deprivation, (d) maximum delay fit residuals with data deprivation

Alternative metrics for quantifying the uniformity of the IPP distribution in the fit domain have been proposed that offer the prospect of additional reduction in the magnitudes of the broadcast GIVEs (Sparks et al., 2003a; Sparks et al., 2003b; Pandya et al., 2007). However, these metrics have not been fully evaluated and have therefore never been implemented.

7 HISTORY OF FIELDED THREAT MODELS

The WAAS undersampled ionospheric irregularity threat model is designed to mitigate threats that arise due to poor sampling. The sampling may be poor either because the spatial distribution of measurements fails to detect the presence of an ionospheric disturbance or because the sampling of an initially well-sampled disturbance deteriorates rapidly over time. When WAAS was commissioned in July 2003, these two facets of the problem were treated separately: a purely spatial threat model table addressed the former, while a purely temporal threat model table addressed the latter (Sparks et al., 2001; Altshuler et al., 2001). From LPV Release 6/7 onward, however, the two facets have been addressed in a single, spatial–temporal threat model. The raw data for each ionospheric threat model have been constructed according to Equation (4). A purely spatial threat model is achieved by setting t_fit,residual to zero; when t_fit,residual > 0, the threat model becomes spatial–temporal.

This section examines the history of the modifications incorporated into each upgrade of the ionospheric threat model. The equations that have governed the construction of each fielded threat model are discussed by Sparks et al. (2026) and summarized in Table A2 of that paper. Table A3 of the current paper summarizes the evolution of the values of the operational parameters used to generate each of these threat models.

7.1 Initial Operating Capability

As discussed by Sparks et al. (2026), the form of σ∼GIVE,ν2, the variance of the distribution that overbounded the vertical delay estimation error, differed slightly in IOC from that given in Equation (1):

σ˜GIVE,v2≡Rirreg−static2(σ˜vmaxcorner)2+max{(Rirreg−staticσdecorrnominal)2,σundersampled,v2}+σROT2,10

where (σ∼vmax corner)2 and R_{irreg–static} were defined, respectively, according to Equations (20) and (17) in the companion paper authored by Sparks et al. (2026), (σdecorrnominal)2 was the delay covariance associated with nearly coincident IPPs, σundersampled,v2 served to bound estimation errors occurring as a result of spatial undersampling, and σROT2 was a variance associated with the maximum rate of change of the level of ionospheric disturbance at the time of the planar fit at the IGP (see Section 7.1.2). Values for σ_{undersampled,v} were supplied by the spatial threat model table.

The tabulation of the IOC ionospheric threat model was based upon observations recorded at 24 of the 25 original WRSs situated throughout the United States and San Juan, Puerto Rico (excluding data from the Honolulu WRS), as depicted in Figure 2. The primary method used to select data for analysis was to scan standard geomagnetic indices to identify days of interest, since the state of the ionosphere is known to correlate strongly with the level of geomagnetic disturbance as indicated by elevated values of the planetary K_p index and depressed values of the D_st index (Datta-Barua et al., 2005). The selected dates initially included January 11, 2000; February 12, 2000; April 6–7, 2000; May 25, 2000; June 8, 2000; July 15–16, 2000; and March 31, 2000. Data recorded on six of these dates were ultimately found to determine the spatial threat model table of the IOC (see Table A1).

Section 7.1.1 describes how values of σ_{undersampled,v} in the spatial threat model table were derived from the tabulation of σ¯undersampled,κ2 as functions of the IPP distribution metrics, discussed in Section 6, characterizing the spatial distribution of fit IPPs around the fit center IGP. Section 7.1.2 discusses the generation of the temporal threat model table used to assign the value of σ_ROT and the Message Type 10 operational system parameters.

7.1.1 IOC Spatial Threat Model

To tabulate the purely spatial threat table of the IOC ionospheric threat model, the estimation of vertical delay and its integrity bound at each IGP assumed that the measurements in each fit were uncorrelated. As discussed by Sparks et al. (2022), this may be considered a limiting case of the more general Kriging approach to delay estimation that was later adopted in 2011 in WFO Release 3A. Fits of slant delay converted to vertical delay were conducted at intervals of Δtfit,W=100 s. To specify the measurement noise covariance (see Equation (7) of Sparks et al. (2026)), IOC WAAS used the definition of the measurement matrix presented by Sparks et al. (2022), with the off-diagonal matrix elements assigned values based on σ_L1–L2,r and σ_L1–L2,s, the standard deviations of the L1–L2 interfrequency bias estimates associated with the measurement receiver and satellite. These standard deviations were set to 1.41 TECU and 1.18 TECU, respectively, and the corresponding variances σL1−L2,r2 and σL1−L2,s2 were assumed to be floor values larger than 99.9% of the actual variances. For the purpose of tabulating the IOC threat model (as noted above), the hysteresis interval during which an irregularity detector could remain in a triggered state following a fit was set to zero. Thus, every threat that failed to trip a local irregularity detector was included in the tabulation.

Figure 7(a) displays the raw data that results from tabulating values of σ¯undersampled,κ2 according to Equation (4) for the IOC data set. Given the relationship between σ¯undersampled,κ2 and |I¯κ−I∼κ|2 defined in Equation (3), it is not surprising how closely the form of Figure 7(a) resembles that of Figure 6(d). Figure 7(b) shows the same data as Figure 7(a) after the overbound has been applied. (Note: Figures 6(c) and 6(d) represent results from a recomputation of the IOC ionospheric threat model performed using the IonoSTAGE software package (Sparks, 2018) developed at the JPL; consequently, the empty pixels in Figure 6(d) [which match the empty pixels in Figure 6(c)] do not exactly match the empty pixels in Figure 7(a) due to minor algorithmic differences and the use of a later version of the set of input observational data.) IOC set the additional conservative inflation factor, K_inflate, in Equation (3) to 1.1. Figure 7(b) includes the impact of this factor.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 7

IOC spatial ionospheric threat model: (a) threat model raw data, (b) threat model overbound

Figures 7(a) and 7(b) use a color bar scale covering the interval [0, 2] to confine attention to the most relevant values of σ_undersampled. Values of σ_undersampled that exceed 2 m give rise to broadcast GIVEs of 15 m or greater, at which point a user’s vertical protection limit generally exceeds the vertical alert limit, making the localizer performance with vertical guidance (LPV) service unavailable. This color bar scale is retained in plots of the raw data and overbounds for all subsequent spatial–temporal threat models.

7.1.2 IOC Temporal Threat Model

The purely temporal threat table of the IOC ionospheric threat model was characterized by a single metric, the time elapsed since the last planar fit; it was constructed by performing fits at 10-s intervals. The planar fit parameters retrieved at a given IGP in a specified epoch were used to evaluate estimates of the vertical delay for all measurements whose IPPs are located within that IGP’s threat domain for some finite duration within the subsequent 800 s. The vertical delay residuals tabulated in the temporal threat model table were defined as follows:

ΔI(t)≡|I¯IPP(t)−I˜IPP(t0)|11

where I¯IPP(t) is the measured slant delay at time t converted to vertical at its IPP in the threat domain of the nearest IGP and I∼IPP(t0) is the estimated value of the vertical delay at the same IPP but evaluated using fit parameters retrieved at the earlier time t₀. As in the spatial threat model, no residual was tabulated when an estimate of vertical delay at the IGP triggered that IGP’s irregularity detector. In contrast to the spatial threat model, however, no data deprivation schemes were applied since it was assumed that the augmentation of error bounds to compensate for poorly sampled irregularities was fully addressed by the spatial threat model.

The results of this tabulation are displayed in Figure 8. As with the spatial threat model table, the delay residual corresponding to each bin was multiplied by an inflation factor K_inflate = 1.1 to provide an additional margin of safety. These data were then used to determine the Message Type 10 operational system parameters that govern the degradation of σ∼GIVE,v over time. The goal of this tuning exercise was to find an overbound of these data that achieved the best performance while still covering the ionospheric threats.

The MOPS (RTCA, 2016) specifies that the adjustment of the broadcast error bounds over time be based on a term (ε_iono in Equation (3) of Sparks et al. (2026)) required to be zero at t = t₀ and increase linearly at least until t = t_message. Since the time dependence of the maximum delay residuals depicted in Figure 8 does not match this requirement, σ_ROT was introduced in Equation (10) to enable the required time dependence to be met.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 8

IOC temporal ionospheric threat model; the blue line labeled ∆I_overbound is an overbound of the delay residuals multiplied by an inflation factor K_inflate ≡ 1.1

Two options defining an overbound were considered: (1) an overbound with a large (initial) constant value and a step-function jump at the initial time of a missed message and (2) an overbound with a smaller initial value but with a ramp that covered all threats at later times. (See Equations (3) and (4) of Sparks et al. (2026), where the degradation provided by these two approaches is specified in terms of the coefficients C_{iono_step} and C_{iono_ramp}, respectively.) Note: no combination of both step and ramp approaches was feasible given the results of the tabulation displayed in Figure 8. Such a combined approach would have been feasible only if the bump in the data that begins at 300 s had occurred more than 100 s later at the initial time of a missed message.

The results of the tuning study indicated that approximately 2% of the time, the ramp approach produced a lower quantized GIVE than the step approach, and coverage improved by 0.2%. Consequently, the ramp approach was selected. Optimal choices for σ_ROT and the Message Type 10 parameters were determined for an overbound expressed as follows:

ΔIoverbound≡KundersampledKinflate[σROT02+Ciono_ramp2(t-t0)2]1/2,12

where K_undersampled = 5.33 as in the inequality expressed in Equation (2). The study concluded that σ_ROT0 and C_{iono_ramp} should be set to 0.30393 m and 0.00075 m/s, respectively (see the ∆I_overbound curve in Figure 8).

7.2 LPV Release 6/7

The first upgrade of the WAAS ionospheric threat model coincided with the fielding of LPV Release 6/7 in September 2007. In Equation (4), the time interval t_fit,residual following each fit epoch over which the maximum magnitudes of fit residuals are tabulated was made finite, eliminating the need to generate distinct spatial and temporal threat models distinguishing between these types of threats. All subsequent ionospheric threat models have been spatial–temporal in nature, addressing both the spatial and temporal types of threats simultaneously.

LPV Release 6/7 introduced nine new WRSs: four in Alaska, two in Canada, and three in Mexico. The IGP mask was augmented to include 131 new IGPs, most of which lay at what was previously the edge of coverage (Sparks et al., 2022). The storm data set used to construct the ionospheric threat model was expanded to include major storms of solar cycle 23 that occurred after the commissioning of WAAS in 2003 (see Table A2). Among these storms were the two that inspired the implementation of the ESD (Sparks et al., 2022): the Halloween storm of October 29–31, 2003, and the extreme storm that followed nearly a month later on November 20–21, 2003. Both the WAAS-only data set and the WAAS+Mexico data set were analyzed in the construction of the threat model, using a much shorter decimation interval than had been the case in IOC: the WAAS-only data were examined at intervals of ∆t_fit,W = 10 s, whereas ∆t_fit,WM, the interval for processing fits of the WAAS+Mexico data, was 30 s, the sampling rate adopted in the tabulation of the Mexican data.

LPV Release 6/7 completely redefined the variance of the overbounding Gaussian distribution used to specify the GIVE at an IGP. The formal error variance for the vertical delay estimate at the IGP for IOC was modified to conform to Equation (6) of Sparks et al. (2026) with c0=(σdecorrnominal)2=(σdecorrtotal)2 . Since the resulting expression for the inflated formal error variance σ∼IGP,V2 (Equation (27) of Sparks et al., (2026)) explicitly includes (Rirreg–dynamicσdecorrnominal)2, it was determined that this σ∼IGP,v eliminated any need for σundersampled,v2 to serve explicitly as a lower bound for the contribution of(Rirreg–dynamicσdecorrnominal)2 to σ∼GIVE,v2.σ∼GIVE,well–sampled,v2 was redefined according to Equation (32) of Sparks et al. (2026) to combine σ∼IGP,v2 with the total maximum antenna bias error term μ_tot,v. The resulting expression for σ∼GIVE,ν2 served to update Equation (22) of Sparks et al. (2026) as follows:

σ˜GIVE,v2=(σ˜IGP,v+μtot,vKHMI−GIVE)2+σundersampled,v2+σROT2,13

where K_{HMI_GIVE} = 5.592 is the constant that defines a confidence interval for a confidence level of (1–2.25×10⁻⁸). The value assigned to K_{HMI_GIVE} matched the GIVE monitor’s fault tree allocation of 2.25×10⁻⁸ per approach as an upper limit on the probability of broadcasting hazardously misleading information (HMI) due to missed detection of an ionospheric disturbance. Since ionospheric irregularities constitute only one of several possible sources of HMI, this probability is required to be less than the 10⁻⁷ value allocated by the fault tree to the entire system (as discussed by Sparks et al. (2026)) to guarantee that a user’s position error is adequately bounded.

In the construction of the LPV Release 6/7 threat model, each table of σundersampled2 values was re-evaluated using the expanded storm data set (see Table A3). The post-fit-epoch interval t_fit,residual over which the evaluation of fit residuals took place was set to a value greater than zero, i.e., the IOC purely spatial threat model was replaced with a spatial–temporal threat model that protected the user not only from spatial threats present in the fit epoch but also from the growth of threats that might occur over the interval t_fit,residual following the fit epoch. The value assigned to t_fit,residual was required to cover the message latency within the user receiver as well as the GIVE computational latency and the system broadcast latency. WAAS CPU loading analyses conducted during IOC, LPV Release 6/7, and LPV Release 8/9.2 determined that IGDs and GIVEs could be computed at three IGPs per second without taxing the system in a critical way. Specifying 330 as an upper bound on the number of IGPs to be included in an IGP mask ensured that the computation of IGDs and GIVEs over the entire IGP mask could be completed within 110 s. Adding this duration to the user receiver message latency of 300 s accounted for setting t_fit,residual to 410 s, a value retained in all subsequent releases to date. In this fashion, the combined spatial–temporal threat model provided protection against the worst-case ionospheric threat that could occur over the intended life of the computed GIVE.

Fits of vertical delay were conducted using the same planar fit methodology that had been used in IOC. As noted previously, however, the measurement noise matrix was set to zero, R_{irreg–dynamic} replaced R_{irreg–static} of IOC in the definition of the irregularity metric at each IGP, t_irreg,hyst and t_irreg,threat were set to 30 and 15 s, respectively, single-station and directional station data deprivation replaced the data deprivation schemes of IOC, and 100-s user smoothing was introduced in the evaluation of fit residuals. A preliminary study of the northern extent of the low-latitude region over which vertical ionospheric delay is no longer well-modeled locally by a plane on days of nominal behavior suggested that the coordinates of the endpoints of the two line segments defining the PPFL, λW,ϕW,λM,ϕM,λE, and ϕ_E, should be assigned the values 26.75°N, 130°W, 18°N, 80°W, 18°N, and 40°W, respectively. The irregularity detector trip threshold was increased to improve system availability; a tuning study of its impact on availability (Pandya et al., 2007) indicated an optimal value of 2.5. The tripping of the ESD eliminated from the threat model severe threats that occurred during the extreme storms of solar cycle 23 (Sparks et al., 2005).

Figure 9(a) shows the raw data that result from tabulating values of σ¯undersampled,κ2 according to Equation (4) using the augmented input data set and the updated processing algorithm. Figure 9(b) shows the overbound of the data in Figure 9(a). When the tabulation of σ¯undersampled,κ2 values was repeated with the WAAS+Mexico data set, the resulting threat model was found to be completely overbound by the WAAS-only component. Thus, no further modification of the spatial–temporal threat model table was required.

Since the overbound depicted in Figure 9(b) included the tabulation of threats that arose in the 410 s following each fit at an IGP, no additional augmentation of σ∼GIVE,v2 was required in LPV Release 6/7 to protect users from the growth of threats over this interval. Consequently, σ_ROT was set essentially to zero. (Actually, it was set to a nominal value of 1 cm and remained in the system as a tunable operational system parameter, available to respond to future needs imposed by the addition of new stations in Mexico and Canada.) The Message Type 10 parameters could now be devoted simply to addressing threats occurring in the period following a missed message.

To protect the user from such threats, the time dependence of σionogrid,v2 continued to be defined by option 2 in Equation (3) of Sparks et al. (2026), but the Message Type 10 parameters in Equation (4) of Sparks et al. (2026) were redefined to accommodate a step-function increase after the 410-s period covered by the spatial–temporal threat model. In contrast to the values adopted in IOC for the Message Type 10 parameters in Equation (4) of Sparks et al. (2026), C_{iono_ramp} was set to zero, and C_{iono_step} was assigned a finite value, chosen to protect users against all ionospheric threats not covered by the 410-s spatial–temporal threat model during the additional 300 s of temporal threat exposure resulting from a missed message. To evaluate this constant, the tabulation period of the spatial–temporal threat model was extended another 300 s to accumulate additional threats over the period of the missed message. Figure 10 displays the regions where the overbound of this threat model exceeded the corresponding regions of the overbound depicted in Figure 9. The worst discrepancy between the two overbounds determined that C_{iono_step} should be set to 0.577.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 9

LPV Release 6/7 spatial–temporal ionospheric threat model: (a) threat model raw data, (b) threat model overbound

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 10

Difference in LPV Release 6/7 between the spatial–temporal ionospheric threat model tabulated over 710 s minus the same tabulated over 410 s

7.3 LPV Release 8/9.2

A relatively minor upgrade of the WAAS ionospheric threat model was fielded a year later in September 2008, in conjunction with LPV Release 8/9.2. This release followed the earlier inclusion in March of four more stations in the WAAS network, two in Canada (Iqaluit & Winnipeg) and two in Mexico (San Jose del Cabo, Tapachula), completing the current WAAS configuration of 38 WRSs. A new IGP mask with 317 IGPs, as plotted in Figure 3(c) of Sparks et al. (2022), was introduced; this mask removed unused IGPs in the Pacific Ocean and south of Mexico at low latitude, while adding new IGPs at high latitude in the northeastern and northwestern regions of WAAS coverage. The same set of storm days upon which the LPV Release 6/7 ionospheric threat model had been based was used to construct the LPV Release 8/9.2 threat model. However, some input files belonging to both the WAAS-only and WAAS+Mexico data set components were first updated to correct the processing error, discussed in Section 2, that was discovered to have impacted the prior computation of supertruth for the storms on July 26–27, 2004; November 8–10, 2004; and January 22, 2005.

The expression for σ∼GIVE,v2 in LPV Release 8/9.2 differed from its predecessor only in that the final term on the right-hand side of Equation (13) was set explicitly to zero. In all subsequent upgrades of the threat model, σ∼GIVE,v2 has remained as follows:

σ˜GIVE,v2=(σ˜IGP,v+μtot,vKHMIGIVE)2+σundersampled,v2.14

Prior to LPV Release 8/9.2, when a fit generated an unphysical, negative estimate for the vertical delay I∼ used to calculate σ¯undersampled,κ2 (see Equation (3)), the value of this estimate was set to zero. In LPV Release 8/9.2 (and in all subsequent updates), this negative fit value was retained in the evaluation σ¯undersampled,κ2, a conservative measure potentially increasing the magnitude of the tabulated values of σ_{undersample,v}.

The most notable algorithmic difference between the LPV Release 8/9.2 threat model and its LPV Release 6/7 counterpart was an adjustment in the position of the PPFL. Moving the PPFL southward decreases the number of measurements recorded at WAAS’s Mexican stations, whose IPPs lie south of the PPFL, or, equivalently, increases the number of measurements available for inclusion in fits of vertical delay at the southernmost IGPs. Thus, this adjustment makes it increasingly likely that the fit domain at each of the southernmost IGPs will circumscribe enough IPPs (> N_min) to conduct a fit. Consequently, an IGP that formerly broadcast a GIVE value of GIVE_MAX (45 m) when the number of fit IPPs was less than N_min can now broadcast a much smaller GIVE based upon such a fit. This, in turn, tends to improve system availability. The disadvantage of moving the PPFL southward, however, is that the additional measurements now incorporated into fits at the southernmost IGPs are increasingly likely to sample regions of the ionosphere at low latitude, where the assumed planar model no longer represents ionospheric behavior accurately. To the extent that these measurements sample ionospheric regions said to exhibit curvature, the values of χirreg2 at southern IGPs may now rise, and the consequent increase in the corresponding GIVEs will tend to impair availability.

Optimizing the PPFL position requires balancing these competing tendencies. Moving the PPFL south of its location in LPV Release 6/7 serves to increase the number of threats included in the ionospheric threat model. It may also increase the frequency of irregularity detector trips at southern IGPs (i.e., χirreg2 is likely to rise when measurements sampling ionospheric curvature are included in the fit). Both responses tend to increase the GIVE values broadcast. Opposing this trend is the reduction in the number of IGPs at low latitude broadcasting GIVE values of GIVE_MAX.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 11

LPV Release 8/9.2 spatial–temporal ionospheric threat model: (a) threat model raw data, (b) threat model overbound

To establish the optimal location of the PPFL, a tuning study (Paredes et al., 2008; Paredes et al., 2009) examined at southern IGPs the dependence of the irregularity metric χirreg2 on the position of the PPFL under nominal ionospheric conditions, both at solar maximum and at solar minimum. This study varied both the value of the irregularity detector trip threshold T_irreg,trip and the coordinates of the points specifying the PPFL location. The study imposed the constraint that the choice of T_irreg,trip and the placement of the PPFL should not substantially degrade performance in the conterminous United States (CONUS) or Alaska at either solar maximum or solar minimum. The study concluded that the PPFL should be moved 5° south of its original position in LPV Release 6/7 but that the irregularity detector trip threshold should remain at 2.5. (Note that in Figures 3(c) and 3(d) of Sparks et al. (2022), the PPFL is erroneously drawn at the location it occupied in LPV Release 6/7 rather than at its new location in LPV Release 8/9.2 and WFO Release 3A, respectively.) Moving the PPFL southward during solar minimum was generally found to provide a substantial improvement in the performance of WAAS throughout Mexico.

Figure 11(a) shows the tabulated raw data for the LPV Release 8/9.2 ionospheric threat model. Figure 11(b) displays the overbound of the data plotted in Figure 11(a). Unlike LPV Release 6/7, the WAAS+Mexico data set in LPV Release 8/9.2 was found to contribute threats that were not covered by the WAAS-only data set. When the period used to tabulate threats was extended another 5 min to accumulate threats over the period of a missed message, the worst discrepancy between the new overbound and that depicted in Figure 11(b) indicated that C_{iono_step} should be reset to 0.457.

7.4 WFO Release 3A

Fielded in October 2011, WFO Release 3 sought to improve system availability by reducing the values of the broadcast GIVEs. To achieve this improvement, WAAS implemented an alternative algorithm for estimating the vertical delay at each IGP (Blanch, 2002; Blanch, 2003; Blanch et al., 2003; Sparks et al., 2010; Sparks et al., 2011a). This algorithm was based on Kriging, a geo-statistical technique, developed originally by the mining industry in the 1950s, that produces a smoothed model of a spatially distributed variable by using a linear least-squares estimator to interpolate values of the variable as sampled by irregularly spaced measurements. Compared to the previous planar fit model, the Kriging delay estimation model generally achieves better agreement with the observed random structure of vertical delay near an IGP by assigning greater weight to the fit measurements whose IPPs lie closest to the IGP.

The introduction of Kriging required an analysis of ionospheric behavior over the full range of possible conditions to determine various tunable Kriging model parameters (as discussed by Sparks et al. (2022)). In addition, it required a reanalysis of the optimal values to be assigned to the trip thresholds and confirmation intervals of the ESD and the irregularity detectors at IGPs. The Kriging parameters to be assigned values were the delay covariance σdecorrnominal associated with nearly coincident IPPs, the delay covariance σdecorrtotal associated with widely separated (uncorrelated) IPPs, and the characteristic decorrelation distance d_decorr (Sparks et al., 2011a). An extensive trade study (Pandya et al., 2012) examined 4 sets of these parameters and 16 values for T_irreg,trip, the irregularity detector trip threshold at each IGP, for a total of 64 system models. The goal of the study was to optimize the availability of the fielded system under moderately disturbed ionospheric conditions while preserving the availability performance achieved using the prior planar fit algorithm under nominal conditions.

Since the operation of the ESD is based on the evolution of the irregularity metric at each IGP and since the implementation of Kriging altered the values of these metrics, it was first necessary to determine, for each set of Kriging parameters, a set of ESD thresholds and confirmation intervals ensuring the ESD would provide a degree of threat protection equivalent to that of LPV Release 8/9.2. The trade study then generated, for each of 64 parameter sets, an undersampled ionospheric irregularity threat model based on the same storm data set used to construct the LPV Release 8/9.2 threat model. The IGP mask differed from that used in LPV Release 8/9.2 only in the inclusion of one additional IGP at (25°N, 55°W). The PPFL was extended 20° eastward, and the WRS at San Juan, Puerto Rico, was now included with the Mexican WRSs assigned to the low-latitude WRS mask governing the removal of IPPs from fits at low-latitude IGPs.

To establish the optimal set of Kriging parameters and irregularity detector trip threshold, the study next simulated system performance for each of the 64 sets of parameters and threat models using a representative set of 14 days of historical system data that spanned the full range of ionospheric behavior over an entire solar cycle. Availability performance was evaluated over a set of four distinct coverage regions: CONUS, Alaska, Mexico, and the entirety of North America. The study concluded that optimal performance was achieved with Kriging parameters σdecorrnominal,σdecorrtotal, and d_decorr set to 0.3 m, 1.0 m, and 8000 km, respectively, the irregularity detector threshold T_irreg,trip set to 3.0, the ESD operational thresholds, namely, the ESD trip threshold T_ESD,trip and the ESD recovery threshold T_ESD,recovery set to 32 and 29, respectively, and the ESD confirmation intervals, namely, the ESD onset confirm interval t_ESD,confirm and the ESD recovery confirm interval t_ESD,recovery, kept at their initial values of 1 h and 8 h, respectively.

Figure 12(a) displays the raw data for the WFO Release 3A ionospheric threat model. Figure 12(b) presents the overbound of the data in Figure 12(a). When threats were tabulated for an additional 5 min to protect the user from the impact of a missed message, the worst discrepancy between the two overbounds required that C_{iono_step} be reset to 0.836. To evaluate in this release the degradation factor that Sparks et al. (2026) define in their Equation (4), option (1) on the right-hand side of their Equation (3) was adopted instead of option (2). Kriging implemented in this fashion was found to improve dramatically both system availability performance and the robustness of the system to challenges posed by ionospheric disturbances, while maintaining system performance under nominal ionospheric conditions.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 12

WFO Release 3A spatial–temporal ionospheric threat model: (a) threat model raw data, (b) threat model overbound

7.5 Release 46-CY16

In the absence of any unusual level of ionospheric disturbance, WAAS generally attains an availability of 100% for the LPV level of service throughout most of CONUS. Prior to Release 46-CY16, however, availability along the coast of California was routinely found to be less than 100%. During the first quarter of 2013, for example, the average availability for the LPV level of service off the coast of California had declined to 99.6%. Even more pronounced was the drop in the availability of the LPV200 level of service a short distance from the coast, where it had fallen as low as 97%. A closer inspection of the data (Sparks & Altshuler, 2014) revealed that temporary loss of availability along the west coast correlated with periodic spikes in the magnitude of the GIVE broadcast by WAAS for the IGP lying in the Pacific Ocean west of Los Angeles at [35°N, 125°W]. The periodicity of these spikes suggested that they were not a repercussion of ionospheric phenomena but rather the consequence of regular variations in the geometry of satellite positions. The time of each spike was discovered to coincide with a time when the fit radius at [35°N, 125°W] exceeded a threshold of 1475 km and simultaneously the RCM exceeded a threshold of 0.575. Such an IPP distribution metric pair identified a region of the threat model where there was a sharp jump in σ_undersampled (see Figure 12(b)).

Fielded in August 2016, Release 46-CY16 sought to improve WAAS availability in coastal regions by introducing into the threat model a dependence on the overall level of ionospheric disturbance throughout the WAAS service volume (Sparks & Altshuler, 2014). This dependence was achieved by introducing the MSD into WAAS operations and constructing the quiet-time branch of the threat model. While the MSD remained in an untripped state, the quiet-time branch provided the σ_{undersampled,v} adjustment to the GIVE at each IGP. Otherwise, this adjustment was provided by the disturbed-time branch, i.e., the same threat model branch depicted in Figure 12(b) for WFO Release 3A.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 13

Release 46-CY16 quiet-time branch of the spatial–temporal ionospheric threat model: (a) threat model raw data, (b) threat model overbound

The disturbed-time branch for this threat model was identical to that shown in Figure 12.

As noted above, the algorithm governing the MSD was similar to that of the ESD, the only difference being lower state detection thresholds and smaller state confirmation intervals. Just as the implementation of the ESD transferred from the threat model to the ESD responsibility for protecting the user from threats that occur during extreme storms, the MSD now triggered the use of the disturbed-time branch to protect the user from threats that occur during moderate storms. Under nominal ionospheric conditions, a much less conservative threat model table could now be used. The quiet-time branch of the threat model was constructed by excluding from tabulation those fit residuals that resulted from fits evaluated when the MSD had tripped. A trade study (Sparks & Altshuler, 2014) concluded that optimal values for the MSD trip threshold T_MSD,trip, MSD onset confirm interval t_MSD,confirm, MSD recovery threshold T_MSD,recovery, and MSD recovery confirm interval t_MSD,recovery were 32, 10 min, 29, and 10 min, respectively.

The same historical storm data set was used to generate the Release 46-CY16 ionospheric threat model as had been used to tabulate its predecessor. Also, the same IGP mask was used. Figure 13(a) displays the raw data for the quiet-time branch of the threat model. Figure 13(b) presents the overbound of the data in Figure 13(a). Note that the quiet-time overbound in Figure 13(b) provides a much less conservative adjustment to the GIVE than does its disturbed-time counterpart (Figure 12(b)), thereby serving to enhance availability whenever the MSD has not tripped. An analysis of WAAS data spanning the interval of January 6, 2011, to November 13, 2013 (Sparks & Altshuler, 2014), concluded that had the MSD been operational with its threshold and confirmation parameters assigned as specified in the preceding paragraph, it would have remained in an untriggered state for over 99.9% of the period. Thus, the implementation of the MSD has promoted a dramatic reduction in the values of the broadcast GIVEs.

7.6 Release 51-CY18

The spatial–temporal ionospheric threat model fielded in September 2018 (Release 51-CY18) was the first to accommodate ionospheric storm data from solar cycle 24 (Sparks & Altshuler, 2016). The ionospheric storms that occurred in solar cycle 24 were uniformly much smaller in magnitude than the largest storms of solar cycle 23 (Sparks et al., 2022). Nevertheless, the storm data of solar cycle 24 included in the Release 51-CY18 threat model computation (storms occurring over the time interval 2011–2015) were found to modify both the quiet-time and disturbed-time branches, due in part to the expansion in the number of WAAS stations. The data upon which all previous ionospheric threat models had been based were collected when the WAAS network consisted of only 25 WRSs. Throughout solar cycle 24, however, this network was comprised of 38 WRSs, including several additional stations located in auroral regions of Alaska and Canada. These new stations provided measurement configurations that differed from those used to generate prior threat models.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 14

Quiet-time branch of the spatial–temporal ionospheric threat model when the tabulation of the Release 46-CY16 threat model is augmented by storm data from solar cycle 24 during 2011–2015: (a) threat model raw data, (b) threat model overbound

Figure 14 shows the impact of solar-cycle-24 data on the quiet-time branch of the ionospheric threat model when data from 2011-2015 recorded at 38 WAAS stations are included in the threat model computation. Comparing Figure 14(b) to Figure 13(b), one finds that the quiet-time branch of the threat model has been significantly degraded, especially in the expanded region covered by red pixels. If this threat model branch had been fielded in Release 51-CY18, system availability would have suffered.

To prevent such a reduction in availability, a new strategy was implemented in Release 51-CY18 to eliminate additional threats from the tabulation of the threat model without diminishing system integrity. The previous threat model algorithm did not account for the fact that, in addition to the protection provided by the various disturbance detectors, a floor value for each GIVE also serves to protect the user from the potentially harmful influence of undersampling. A floor value is broadcast whenever the computed value of the GIVE at an IGP falls below the floor value specified for that IGP. The floor value at any IGP has never been set in WAAS to a value less than 3 m. The threat model algorithm in Release 51-CY18 was upgraded to exclude threats using not only the ESD, MSD, and IGP irregularity detectors but also to exclude a threat when it is covered by the user ionospheric vertical error (UIVE) at the threat IPP, as determined by interpolating the GIVE floors at the corners of the cell in which this IPP is located. In other words, to tabulate a threat, σ¯undersampled,κ2 had to satisfy both of the following conditions (see Equation (3)):

σ¯undersampled,κ2>015

σ¯undersampled,κ2+σ˜κ2>σUIVE−floor,κ2,,16

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 15

Release 51-CY18 spatial–temporal ionospheric threat model: (a) quiet-time raw data, (b) quiet-time overbound, (c) disturbed-time raw data, (d) disturbed-time overbound

where σ∼κ2 is the inflated variance of the delay estimate at that IPP and σUIVE_floor,κ2 is the variance of the UIVE floor at the user’s (threat) IPP, calculated by interpolating the GIVE floor variances at the neighboring IGPs. (Recall that these variances depend upon the magnitude of the quantized GIVE as specified in the MOPS (RTCA, 2016).) This method of restricting the threats tabulated in a threat model has been designated UIVE floor threat culling.

Figure 15 displays the raw data and overbounds for both branches of the Release 51-CY18 ionospheric threat model. In addition to excluding low-latitude IGPs that had been present in the tabulation of the WFO Release 3A threat model, the IGP mask for this threat model modified the GIVE floor values at [15°N, 100°W] and [15°N, 105°W] to change from 15 m to 3 m.

To produce overbounds of the raw threat data tabulated over 710 s and 410 s, the same procedure was initially followed as in WFO Release 3A. When the difference between these tables was calculated, however, the worst discrepancy defining C_{iono_step} was calculated to be 1.1213 m, which exceeded the range for C_{iono_step} allowed by the MOPS (RTCA, 2016). (The previous 0.836 value of C_{iono_step} was retained.) To preserve system integrity, the disturbed-time branch tabulated over 710 s (displayed in Figure 15(d)) was simply substituted for the corresponding branch tabulated over 410 s. This highly conservative decision did not affect the tripping of the MSD and had little impact on availability in solar cycle 24, since the disturbed-time branch of the ionospheric threat model was seldom accessed during this period. Thus far in solar cycle 25, however, MSD trips have been far more frequent, resulting in a larger impact on availability. The decision to substitute the disturbed-time branch tabulated over 710 s for the same tabulated over 410 s needs to be re-evaluated.

7.7 Release 62-CY23

In May 2023, a minor upgrade of the WAAS ionospheric threat model was fielded in Release 62-CY23. The quiet-time and disturbed-time branches of this upgrade differ only slightly from those depicted in Figure 15. This upgrade was motivated primarily by the discovery of WAAS observations, recorded near the southern border of the WAAS coverage region on June 1, 2013, and June 29, 2013, that did not pass the UIVE bounding (UB) analysis required to justify the safety and integrity of WAAS. The UB tool used to perform this analysis calculates residual differences between vertical delay estimated at a given set of test IPPs (as determined by system observations) and the corresponding set of actual slant delay measurements converted to equivalent vertical delay at the IPPs. The estimated vertical delay value and the UIVE evaluated at each test IPP location are determined by interpolating, respectively, the IGDs and the GIVEs, at the IGPs nearest the test IPP, according to the MOPS bilinear interpolation algorithm (RTCA, 2016). In this analysis, these IGDs and GIVEs are determined from fits that exclude the test point in question. Furthermore, the values of σ∼GIVE,v used to define the GIVEs for this analysis are evaluated according to Equation (1) with the threat model contribution set to zero. The UIVEs are used to normalize the residuals, and the distributions of the resulting normalized residuals are tabulated in histograms on which various bounding analyses are performed. The results of this study determine whether, under both nominal and storm conditions, the GIVE algorithm generates UIVEs sufficient to bound the user’s ionospheric error.

Figure 16 shows the type of event that caused the UB analysis to fail when analyzing observational data recorded in June 2013. Figure 16(a) displays tracks of slant delay, converted to equivalent vertical delay according to the instantaneous value of the satellite elevation angle, that have been derived from observations of GPS signals emitted on June 1, 2013 by satellite PRN 5 (note: each GPS satellite is identified by its pseudorandom noise (PRN) number), and recorded by the three receivers at the Coyuca station in Mexico, roughly 20 miles northwest of Acupulco. Tracks are displayed as observed by the system in real time (threads A and B only) and as determined by post-processing that has honed the data (all three threads). Note the precipitous decline in the delay that occurs near 3:00 Coordinated Universal Time (UTC). This is a manifestation of a deep and narrow TEC trough that occurred off the coast of Mexico and generated an ionospheric threat that was not adequately addressed by any of the system’s means of threat mitigation.

Figure 16(b) shows the values of the equivalent vertical delay for all WAAS measurements recorded in the epoch at 3:10 UTC. The only measurement that detects the presence of the TEC trough is represented by the turquoise marker located in the red box, identifying the threat domain of the IGP located at the center of the box. Note, however, that the turquoise marker lies below the PPFL, and this measurement is therefore excluded from the fit of equivalent vertical delay centered on the IGP. As a consequence, the system was, at this time, ignorant of the presence of the TEC trough near Acapulco.

Various possible modifications to the existing WAAS means of threat mitigation were proposed to address this problem, such as tuning the parameters governing Kriging, tuning the tripping of the local irregularity detector, and tuning the location of the PPFL. The simplest solution, however, was simply to increase the values of the GIVE floors assigned to IGPs at the southern edge of the WAAS service volume. In Release 51 – CY18, GIVE floors had been assigned values uniformly of 3 m throughout the coverage region. Since ionospheric threats tend to be distributed as a function of geomagnetic latitude (see Section 7 of Sparks et al. (2026)), it was decided to keep the GIVE floors constant within geomagnetic latitude bands of 5° width (depicted in Figure 17), as determined by the methodology of Bowring (1976).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 16

The TEC trough that occurred off the coast of Mexico on June 1, 2013: (a) measured slant TEC tracks of TEC, converted to vertical delay, for three GPS signals propagating from PRN 05 to the Coyuca station, as recorded by the system in real time and post-processed, (b) distribution of IPPs showing the equivalent vertical delay for slant delay measurements at 3:10 on June 1, 2013

In Figure 16(b), the red box indicates the threat domain of the IGP nearest the IPP of the only measurement to sample the TEC trough. The two contiguous purple line segments define the position of the PPFL.

To determine the magnitude of the necessary increase in the values of the GIVE floor at IGPs in southern Mexico, a study was conducted to assess the minimum value that each GIVE floor could take on, designated the critical GIVE floor, to mitigate all threats that arise within that IGP’s threat domain for the observational data set used to tabulate the threat model. Setting the GIVE at a given IGP to the critical GIVE floor at that IGP eliminates any need to augment the GIVE using an ionospheric threat model. Critical GIVE floors were computed by two distinct methods. The initial method took into consideration the interpolation required to evaluate the UIVE: whenever a threat whose IPP lay in a given grid cell was not adequately bounded by the UIVE, GIVE floors at the four corners of the cell were raised sequentially until the threat was bounded. The second method made use of the fact that each threat incorporated into the threat model is determined by an estimate of the vertical delay at the threat IPP based on only a fit centered on the nearest IGP (no interpolation is involved). The critical GIVE floor can then be determined simply by raising the GIVE floor at the fit IGP until its magnitude is sufficient to cover the threat.

When data covering the June 2013 events in question were included in the analysis, the results revealed that the IGPs in the southernmost bands needed to have their GIVE floors raised to at least 4.5 m. Initially, the GIVE floors in the first two geomagnetic bands (bands 20°–25° and 25°–30° in Figure 17) were assigned values of 4.5 m and 3.6 m, respectively. With this configuration, however, the UB analysis continued to fail. Upgrading these values to 6.0 m and 4.5 m, respectively, enabled the UB analysis to pass.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 17

IGPs in 5° bands of geomagnetic latitude, where each band is assigned a constant GIVE floor

In addition to altering the GIVE floor values assigned to the southernmost IGPs, the ionospheric threat model of Release 62-CY23 differs from its predecessor in another important respect – the storm set used to generate it. For Release 62-CY23, an expanded set of criteria has been used to determine ionospheric storm days of interest. The motivation for this expanded set was both to tighten the criteria identifying days of interest and to accommodate situations in which the previous criteria missed ionospheric disturbances of potential interest. In addition to the standard geomagnetic disturbance indices K_p and D_st used previously to identify days when the ionosphere exhibits a high level of disturbance, WAAS now takes into consideration IGP irregularity detector trips and the magnitudes of drops in coverage. These new criteria were first examined in the selection of days of interest to be included in the tabulation of the Release 51-CY18 threat model. These criteria were evaluated as a possible means of identifying days of interest when the geomagnetic indices do not manifest a high level of disturbance. Ultimately, however, all the storms selected for that threat model had values of K_p that attained at least a magnitude of 7-.

To establish a disturbance category index (DCI), the IGPs comprising the WAAS grid are grouped into six regions that correspond roughly to CONUS, Alaska, Canada, Mexico, the Atlantic Ocean, and the Pacific Ocean (see Figure 18). On a given day, the level of ionospheric disturbance in each region is characterized independently by the magnitude of any drops in coverage and a statistical metric ξ representing the fraction of the day during which the region contained one or more IGPs whose irregularity metric exceeded 1.0 (a lower threshold than the detector trip threshold of 3.0). The definition of the DCI makes use of prior calculations that have determined a bound for ξ in each region, designated ξ_1%, computed such that precisely 1% of the days in a representative data set exhibit values of ξ that exceed ξ_1%. The data set chosen to determine these bounds was the set of measurements recorded at WAAS stations over the two-year duration 2011–2012, a period of generally modest ionospheric activity. The DCI also assigns a low but non-negligible priority to days of interest when the LPV and LPV200 levels of service exhibit significant drops in coverage in response to ionospheric activity characterized by K_p ≥ 4.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 18

Geographic regions used to define the DCI in Release 62-CY23 identifying days of interest to be included in the construction of the WAAS ionospheric threat model

Table A2 identifies the attributes associated with each value of the DCI that can be assigned to a storm day. Based on five flags, this index gives the highest priority to days on which the maximum K_p index is greater than or equal to 7- and the minimum of the D_st index drops below –100 nT. The index gives secondary priority to days on which Nξ,mid-low>1 and Nξ,high>0, where Nξ,mid-low and Nξ,high are respectively, the number of regions at middle and low latitude (CONUS, Mexico, Atlantic Ocean, Pacific Ocean) and the number of regions at high latitude (Alaska and Canada) in which the ξ value for a region exceeds the ξ_1% assigned to that region. The DCI also assigns a lower but non-negligible priority to significant drops in coverage for the LPV and LPV200 levels of service. Three levels of coverage (<100% <99%, <95%) are monitored, both for each level of service in CONUS and for each level of service in Alaska. Drops in coverage are flagged when N_{coverage,CONUS}, the number of coverage-drop conditions satisfied in CONUS, is greater than 1, and when N_{coverage,Alaska,} the number of coverage-drop conditions satisfied in Alaska, is greater than 1 (these coverage-drop conditions are evaluated only for days when K_p becomes greater than or equal to 4).

The implementation of the DCI provides a means of ranking each day throughout solar cycle 24 and beyond according to its level of ionospheric disturbance, thereby providing a means to determine the days whose contributions to the computation of an ionospheric threat model are of potential interest. (Not all of the required data remain available to enable evaluation of DCIs for solar-cycle-23 storm days.) In an analysis distinct from that used to determine ξ_1% for each geographic region, 137 days were studied over the time interval 2011–2018 to develop criteria for selecting days of interest. Days with DCIs of less than 12 were found to have a low impact on the threat model. Therefore, a value of 12 has been selected as the minimum DCI for identifying an ionospheric day of interest whose measurement data are to be considered for inclusion in a threat model computation. Note that the DCI has been designed such that storms whose K_p or D_st values transcend respective assigned thresholds are automatically selected for inclusion in the computation.

To determine the list of storm days to be analyzed for the construction of the Release 62-CY23 threat model, a single-day threat model was generated for each day in 2011–2018 when the DCI was greater than or equal to 12. Any single-day threat model with critical points whose σ_{undersampled,v} values rose above the corresponding points of the previously fielded threat model (Release 51-CY18) was added to the list. Days of interest with critical point σ_undersampled values just below those of the fielded threat model were also included in the list, to allow for addressing the possibility that future algorithm changes might result in such days adding critical points to the threat model. Specifically, days of interest with the following pixel differences were included in the list of storm days associated with the threat model: absolute pixel differences <0.35 m and relative pixel differences <25%. For Release 62-CY23, this methodology added four new storm days (June 29, 2013; February 27–28, 2014; December 20–21, 2015; March 6–7, 2016) to the list of storm days processed in the construction of the Release 51-CY18 threat model.

8 SUMMARY

This paper has sought to preserve a description of the rationale and the decisions that have governed the construction and evolution of the WAAS undersampled ionospheric irregularity threat model from its origin in IOC to the present. The threat model is derived from observations of GPS signals recorded by WAAS receivers, processed to produce ionospheric truth data designated supertruth. To ensure that each ionospheric threat model adequately protects users at low latitude, the observational data supplied by WAAS stations have been augmented by historical data recorded at Mexican sites during solar cycle 23, a period prior to the installation of WAAS receivers in Mexico. The processing of observational truth data to generate a threat model takes into account the threat to position estimate accuracy posed by the smoothing of delay measurements performed by WAAS user receivers. Since the actual observations used to construct a threat model do not necessarily capture the worst possible placement of ionospheric threats relative to the sampling measurements, the observational data used to construct a threat model are augmented by simulated configurations of measurements derived by systematically removing various subsets of observations from sets of actual observation, a technique designated data deprivation. In the tables constituting the threat model, the two metrics used to characterize a distribution of measurement IPPs incorporated into each fit of vertical delay are the radius of the fit domain and the ratio of the fit IPP centroid radius to the fit radius.

This paper has documented the changes introduced in each upgrade of the ionospheric threat model from WAAS’s commissioning on July 10, 2003, to the threat model fielded most recently in 2023.

HOW TO CITE THIS ARTICLE:

Sparks, L., Altshuler, E., Blanch, J., Walter, T., McCord, E., & Griffin Sanchez, R. (2026). WAAS and the ionosphere – a historical perspective: threat model evolution. NAVIGATION, 73. https://doi.org/10.33012/navi.758