Ionospheric spatial decorrelation assessment for GBAS daytime operations in Brazil

Extensive ionospheric studies were conducted to support the initial phase of system design approval for the existing SLS-4000 GBAS installed at Antonio Carlos Jobim International Airport (formerly Galeão International Airport) (GIG) in Rio de Janeiro, Brazil. This paper focuses on determining the broadcast value of the standard deviation of vertical ionospheric gradients (or 𝜎 vig ) that is required to bound ionospheric spatial gradients in Brazil under nominal conditions during daytime hours. The time-step method is useful for gaining sufficient samples at distances less than the physical separation distance of ground stations and was utilized to estimate ionospheric spatial gradients. A new method called “geometric similarity” was developed to estimate ionospheric temporal gradients and evaluate the temporal effect added to the bounding 𝜎 vig values. As a result, a 𝜎 vig of 13 mm/km, including a temporal gradient contribution of approximately 2 mm/km, is conservative enough to bound ionospheric spatial decorrelation for daytime GBAS operations in Brazil.


INTRODUCTION
The Ground Based Augmentation System (GBAS) is designed to enhance the performance of Global Navigation Satellite Systems (GNSS) to enable precision approach and landings for civil aviation. A ground facility for GBAS is typically equipped with four GNSS reference receivers at known survey positions within the property of a particular airport. Based on the difference between actual and GNSS-calculated positions, the ground facility produces differential corrections as As part of aviation research and development in Brazil, the SLS-4000 was installed at Antonio Carlos Jobim International Airport (formerly Galeão International Airport) (GIG) in Rio de Janeiro. Ionospheric activity in equatorial regions (low-latitude regions within 25 degrees latitude of the geomagnetic equator), including Brazil, is known to be significantly more variable and intense than that encountered (and extensively studied) in midlatitude regions (Abdu et al., 2003;Venkatesh, Tulasi Ram, Fagundes, Seemala, & Batista, 2017). However, the SLS-4000 originally was configured with ionospheric parameters and threat models for the conterminous United States (CONUS), which is almost all in the mid-latitudes. Thus, the CONUS threat model and ionospheric parameters needed to be re-examined before they could be applied to GBAS in Brazil.
The Brazilian GBAS ionospheric threat model study was carried out by an international and interagency team composed of experts on the ionosphere and GBAS starting in October 2013 (Yoon et al., 2017). This study utilized dualfrequency GNSS data collected on 123 active ionospheric days during the peak of Solar Cycle 24 (April 2011-March 2014. Over 1,000 anomalous gradients were identified from post-processed data and were fully validated as being due to ionospheric behavior. Among them, 59 validated gradients exceed the maximum of 412 mm/km observed in CONUS. Of these, 31 gradients were observed that exceeded 500 mm/km, with a maximum of 850.7 mm/km. However, those gradients which exceeded the CONUS threat model bounds occurred during local nighttime, when disturbed ionosphere and equatorial plasma bubbles are known to occur. The few events that occurred outside local nighttime did not exceed the CONUS threat model bounds. The results of Yoon et al.'s (2017) study suggest that, while additional mitigation would be needed for Category I (CAT I) GBAS to be provided at nighttime, the existing GBAS with the CONUS threat model might be sufficient to support operations during daytime. As an extension of the ionospheric threat model study (Yoon et al., 2017) that was based on Brazilian ionospheric data from 2011 to 2014, additional days of severe ionospheric storms during March 2015 and September 2017 and well-known storm days during the previous Solar Cycle 23 in 2003 and 2004 were examined (Lee & Pullen, 2018). This study also confirmed that the CONUS threat model bounds all gradients that have been observed during daytime hours.
While no additional ionospheric threat mitigation is needed beyond the existing GBAS during daytime hours, the nominal standard deviation of vertical ionospheric gradients ( vig ) over Brazil should be assessed to make sure that it bounds ordinary low-latitude ionospheric conditions. Small user errors remain after GBAS corrections due to ionospheric spatial decorrelation between reference station and users. This spatial decorrelation error is taken into account when the user computes protection levels (PL) which bound the true position error to the required integrity probability. This is done in the calculation of PL by applying the broadcast vig , which is selected to meet the integrity requirements of CAT I precision approaches. The vertical protection level (VPL) for a nominal or faultfree case is given in a form by where is a multiplier determined by the probability of a fault-free missed detection, , is the vertical position component of the weighted-least-squares projection matrix for satellite , and is the standard deviation of a normal distribution that overbounds the true range domain error distribution under the fault-free hypothesis (RTCA, 2008). The broadcast vig is used in the computation of .
In mid-latitude CAT I GBAS operations, the value of vig that bounds nominal ionospheric spatial gradients has been determined to be 4.0 mm/km in CONUS (Lee, Pullen, Datta-Barua, & Enge, 2007) and 2.07 mm/km in Germany (Mayer, Belabbas, Jakowski, Meurer & Dunkel, 2009). The actual value that is broadcast by the SLS-4000 is the Root-Sum-Square (RSS) of 4.0 mm/km and a separate value of 5.0 mm/km that was conservatively chosen to bound anomalous tropospheric gradients (van Graas & Zhu, 2011). Therefore, the minimum broadcast value of vig in CONUS is (4.0 2 + 5.0 2 ) 0.5 ≅ 6.4 mm/km. The same approach can be taken in Brazil, but the larger spatial gradients created by equatorial ionospheric behavior will require a higher value of the ionospheric component of vig . As an example, the value of vig in Japan located in low-and mid-latitude regions was estimated to be 6.02 mm/km (Yoshihara, Saito, & Fujii, 2010).
Prior work estimated ionospheric spatial gradients using two existing methods, the station-pair method and the time-step method (Lee et al., 2007). Because of its architectural resemblance to the GBAS ground station and aircraft configuration, the station-pair method is a more intuitive means to measure ionospheric spatial gradients. However, the sparsity of GNSS reference stations in Brazil makes it difficult to use this method to estimate ionospheric gradients over short distances because the distance over which the spatial gradients are estimated depends on the physical separation of stations. Unlike the station-pair method, the time-step method groups a single satellite and a single receiver as a pair and obtains the spatial separation of interest by adjusting the time window of observation.
Thus, we utilize the time-step method to estimate vig in Brazil. However, the time-step method introduces temporal errors, meaning that the estimated gradients would be a mixture of spatial and temporal variations. In this study, we also investigate the effect of temporal decorrelation errors added to the estimated vig by extracting temporal gradients from the total gradient estimates.
This paper determines the minimum value of vig which is required to bound ionospheric spatial gradients at Galeão airport under nominal conditions during daytime hours. Section 2 describes the datasets used to estimate nominal ionospheric spatial gradients in Brazil. Section 3 explains the time-step method that was used to estimate ionospheric spatial decorrelation. Section 4 gives estimates of the bounding nominal vig over Brazil during early morning, daytime, and early evening hours and recommends a specific value of vig that is valid over a recommended local daytime window. Section 5 explains how to translate this general guidance regarding the duration of the daytime window into specific times to begin and end GBAS operations that can be derived (as a function of season) from GPS time. In Section 6, the temporal gradient effects included within the bounding vig are assessed by estimating temporal gradients via a newly proposed "geometric similarity" method. Finally, Section 7 draws conclusions along with recommendations for future work.

Datasets from the Brazilian GNSS networks
As described above, the broadcast value of vig is required to bound nominal ionospheric gradients, meaning those that may be present under all but anomalous (plasma bubble or storm-driven) conditions. Thus, "nominal" includes both ionospherically quiet and active conditions, because GBAS cannot distinguish between these conditions in real time (Lee et al., 2007). Ionospheric gradients under the most severe conditions, defined as anomalous, are included in a GBAS threat model and are usually (but not always) detected by GBAS monitors (Lee, Seo, Park, Pullen, & Enge, 2011;Pullen, Park, & Enge, 2009). While the goal of this work is to determine the value of vig during daytime, it utilizes data collected for the Phase I Brazilian threat model study (Yoon et al., 2017) that observed extremely large gradients generated by equatorial plasma bubbles at nighttime. Several studies discovered that ionospheric irregularities initially observed on a previous night or near sunrise time can persist during daytime hours due to the long lifetime of irregularities (Fukao, (Yoon et al., 2017) [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] 2010). The results of these studies suggest that these nighttime ionospheric events could affect daytime ionospheric conditions. In order to include all nominal and ionospherically active days (that are not anomalous during the daytime) in daytime vig estimation, we selected dates from the Phase I study, which includes the most disturbed ionospheric conditions at nighttime and thus is very likely to bound all ionospheric conditions in daytime.
The Phase I study in Yoon et al.'s (2017) paper analyzed data from several networks of ground stations in South America, including the Brazilian Network for Continuous GPS Monitoring (RBMC), the Low-Latitude Ionospheric Sensor Network (LISN), the Concept for Ionospheric scintillation mitiGAtion for professional GNSS in Latin America (CIGALA), the ICEA Septentrio and Trimble receiver network, the Integrated Positioning System for Geodynamic Studies (SIPEG), and the Honeywell SLS-4000 receivers installed at and near Galeão airport (GIG). Figure 1 shows the locations of the approximately 185 GNSS reference stations used in this study (as of 2014), color-and symbol-coded by the network that each station belongs to. The densest region of reference stations corresponds to the most populated regions of Southeastern and Coastal Brazil.
The total of 123 days of data from April 2011 to March 2014 in the peak of Solar Cycle 24 were selected to search for the worst ionospheric gradients and determine the upper bound of the threat model in Yoon et al.'s (2017) study. This selection was done by utilizing pre-existing indicators of space weather intensity, including planetary K (K p ) and disturbance/storm time (D st ) values. While these geomagnetic indices characterize global ionospheric disturbances fairly well, they fail to completely capture the local ionospheric distortions for low latitudes. Thus, the L1-band amplitude scintillation index S 4 was also utilized to identify post-sunset equatorial ionospheric instabilities.

Date selection
Strong ionospheric scintillation is a feature of low-latitude ionospheric behavior that distinguishes it from midlatitude behavior, where scintillation is much rarer. Scintillation events are known to occur mostly at night, after local sunset, but, as mentioned previously, they can extend toward the early morning and late-afternoon hours. Thus, the days for the analysis were selected by a combination of high levels of daytime or nighttime scintillation (as measured by the Scintillation Occurrence Rate or SOR) and/or severe values of the D st index. SOR is the ratio of the number of points in a time period in which the S 4 amplitude scintillation index exceeds a threshold to the total number of points in that time interval and is defined as follows (Guo, Liu, Zhao, & Wang, 2017): where N[all] is the total number of points observed within a certain time period, and N[S 4 ≥ 0.2] is the number of points in that time interval in which the S4 amplitude scintillation index exceeds 0.2, indicating a significant degree of amplitude scintillation. From the database of 123 days identified in the Phase I study, two ionospheric datasets were selected based on the highest daytime SOR and the highest nighttime SOR, respectively. The lists of the selected days, their data categories from the Phase I study, SOR with rankings, and bounding vig results are shown in Table 1 (the top five daytime SORs) and Table 2 (the top six nighttime SORs).

IONOSPHERIC GRADIENT ESTIMATION METHOD
Two methods, the station-pair method and the time-step method, were developed and used to estimate anomalous ionospheric spatial gradients in previous work (Datta-Barua et al., 2010;Lee et al., 2007). The station-pair method is illustrated in Figure 2 (left). Estimates of ionospheric slant delay are made by two ground stations (in known locations) from the same GPS or GNSS satellite at the same time, and the spatial gradient is calculated simply as the difference between these two delays divided by the baseline distance between the two stations. Note that distances within the presumed ionospheric thin shell are deliberately not used, as the thin-shell height is much more uncertain under anomalous conditions. For a similar reason, slant ionospheric gradients are used instead of applying the ionospheric obliquity factor (derived from the thin-shell model at a specific height) to calculate elevation-independent vertical or "zenith" gradients.
The time-step method is also used to estimate ionospheric gradients, and it is shown in Figure 2 (right). It only requires one ground station, and it estimates gradients over the path in space created by satellite motion from ionospheric pierce point (IPP) at 1 to IPP at 2 over a F I G U R E 2 Station-pair method for estimating anomalous ionospheric gradients (left) and time-step method for estimating nominal ionospheric gradients and validating anomalous ionospheric gradients (right) [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] short time interval Δ = 2 − 1 . Here, the distance over which the change in ionospheric delay is divided by is based on the presumed IPP separation distance based on the ionosphere being concentrated in a thin shell at a particular altitude (typically 350 km). As noted above, this is a poor assumption under anomalous conditions, but it is more accurate under nominal conditions. Also, because gradients are estimated over a time interval instead of at a single time, the results include temporal as well as spatial changes in ionospheric delay, meaning that they overestimate or underestimate the spatial gradients that are the desired results. Given these limitations, the time-step method is used only for estimating gradients on nominal days (as done in Lee et al.'s 2007 paper, accepting that the results will be conservative), validating anomalous gradients estimated using the station-pair method (Yoon, Kim, & Lee, 2017), and as a backup method for estimating anomalous gradients when the density of ground stations is insufficient to provide good estimates from the stationpair method.
In this study, the time-step method was used to estimate ionospheric gradients, since the density of South American ground stations from which data was collected and analyzed is insufficient to provide good estimates from the station-pair method. This method is useful for gaining sufficient samples at distances less than the physical separation distance of ground stations, especially when the distances between nearby stations are more than 100 to 200 km apart. The time window (Δ = 2 − 1 ) of one minute was used in this study to obtain the spatial distance of interest (i.e., short distances of less than 25 km representing the effective separation between aircraft approaching an airport to land and the GBAS ground facility at that airport). However, ionospheric temporal decorrelation over one-minute intervals is included in the spatial gradient estimates, and this will be further examined in Section 6.  Figure 3 shows the spatial decorrelation result for the day with the highest daytime SOR value (March 7, 2012, which is ID 1 from Dataset 1 in Table 1 differential delays into bins, and the color scale of each pixel indicates the number of measurements counted in each (x, y) bin. The absolute value of the vertical delay difference increases fairly linearly as the IPP separation distance increases. From these results, the differential delays were divided by the corresponding IPP separation distances to obtain vertical ionospheric gradients. The Complementary Cumulative Distribution Function (CCDF, representing 1 minus the CDF function and thus showing the probability of exceeding a given x-axis value) of normalized vertical ionospheric gradients is shown in Figure 4 on a semilogarithmic scale. The vertical gradients are normalized by removing their mean ( vig ) and dividing them by their standard deviations ( vig ) derived from the observations shown in Figure 3. The CCDF of the normalized actual gradient distribution and that of the standard normal (Gaussian) distribution are denoted by the red solid curve and the black dashed curve, respectively. Figure 4 shows that the actual distribution derived from the observations is not bounded by the CCDF of the standard normal distribution between probabilities of approximately 7 × 10 -2 and 6 × 10 -5 and thus has non-Gaussian tails. Because GBAS users assume a zero-mean Gaussian distribution of errors in the computation of protection levels, the nominal sigma ( vig ) of a zero-mean Gaussian distribution must be inflated to cover the non-Gaussian tails of the actual distribution. Here, we have inflated the nominal sigma until the CCDF of normalized vertical gradients does not exceed that of the standard normal distribution. The infla-F I G U R E 5 σ vig overbound results from the time-step method and Brazilian GNSS network data on the highest daytime SOR day (March 7, 2012) [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] tion factor needed for the March 7, 2012, data is 1.125. The CCDF of vertical gradients re-normalized by the inflated sigma (1.125 × vig ) is represented by the blue dash-dotted curve in Figure 4, and it mostly falls well below the CCDF of the standard normal distribution. This indicates that the inflated Gaussian distribution bounds the empirical distribution and has with significant margin for most CCDF probabilities. However, any smaller inflation factor will not bound the standard normal CCDF where they come together around a probability of 4 × 10 -3 ; thus 1.125 represents the minimum inflation factor for this data set.

overbound estimation method
In this study, the " vig overbound" estimation method presented in Lee et al. (2007) was used to assess the nominal vig value for Brazil. The mean and standard deviation of vertical ionospheric gradients in each bin of IPP separation distance are computed, interpolated to the distances corresponding to each gradient, and are then used to normalize the gradients. Based on the distribution of normalized ionospheric gradients, an inflation factor ( ) of 1.125 is then determined empirically to achieve the bounded condition shown in Figure 4. Lastly, the overbounding vig is computed as | vig | + * vig for each bin. Figure 5 shows the vig overbounding result on March 7, 2012. The estimated vig overbounds (the curve with asterisks) and the one-sigma values (the curve with circles) are below 10 mm/km and 8 mm/km, respectively, at separation distances less than 6 km. While the estimates at distances less than 4 km are not larger than 10 mm/km for this day, those at the shorter IPP separation distances are not reliable because of the relatively large effects of temporal decorrelation and residual biases (due F I G U R E 6 Summary of σvig overbound results for each Data ID in Datasets 1 and 2 [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] to carrier-phase leveling error) at these distances. These effects degrade the results more when divided by shorter IPP distances. Estimates at distances larger than 6 km are not reliable enough to be used either because they were obtained using an insufficient number of samples for reliable statistics (only about 13% of the total data points are for IPP distances larger than 6 km). This is supported by observing the increasing vig estimates as a function of IPP distance, while the opposite trend is expected due to the reduced impact of temporal correlation and residual biases as shown in Lee et al. (2007). Nevertheless, the results obtained at IPP separation distances of 4-6 km, which correspond to the range of distances between the location of the precision approach decision height and the GBAS ground facility, support the conclusion that vig is conservatively bounded by 10 mm/km on March 7, 2012.
Overbounding vig values were estimated for all days listed in Table 1 and Table 2 using the same process conducted for the case of March 7, 2012. These results are shown in the following subsection.

overbound results for selected datasets
The analysis in this study examined the two ionospheric data sets from the Phase I threat analysis database explained in Section 2 using an approach similar to that of Lee et al.'s (2007) study but using the time-step method to derive spatial gradient estimates in place of the stationpair method, which requires a denser network of ground receivers. A summary of the vig overbound results for each ID from Datasets 1 and 2 is shown in Figure 6. The data sets selected based on the daytime SOR and nighttime SOR produced bounding vig values in the range of 8-10 mm/km and 10-13 mm/km, respectively, when daytime was defined to exist between 6 a.m. and 6 p.m. local time. The " vig " value here refers only to the nominal ionospheric component and not to the RSSed tropospheric component included within GBAS.
One of the days identified, October 1, 2012, in Dataset 2 (ID 6) in Table 2, turned out to be an anomalous day because it included gradients exceeding 100 mm/km (50-100 times larger than a typical one-sigma value of vig ). Figure 7 shows the slant ionospheric delays (top) and vertical ionospheric gradients (bottom) observed from four nearby station-satellite pairs between 5 a.m. and 7 a.m. local time on that day. Gradients larger than 100 mm/km were observed not only right after the sunrise time (5:30 a.m. local time) from the PPTE-PRN6 (downward triangle), PRMA-PRN22 (circles), and PRU2-PRN6 (squares) stationsatellite pairs, but also after 6 a.m. local time from the PISR-PRN22 (asterisk) station-satellite pair. Similar patterns of slant ionospheric delays experienced by multiple station-satellite pairs demonstrate that the observed anomalous gradients are real. Thus, this day was considered as anomalous and was excluded from the days used to estimate the nominal vig value. (zenith) differences in ionospheric delay as a function of IPP separation distance for data sets starting at 6 a.m. local time (marked with blue circles) and 5 a.m. local time (marked with red diamonds) on this day. Much larger variations appear in the data that starts at 5 a.m. local time. This was caused by an ionospheric disturbance occurring prior to sunrise (at 5:50 a.m. local time) and affecting multiple satellites viewed from station POAL (Porto Alegre). Figure 9 shows the estimates of slant ionospheric delays (top) and vertical ionospheric gradients (bottom) observed from station POAL viewing PRN3 (downward triangle), PRN6 (circle), and PRN22 (square). The rapid variations in slant delays generated large gradients, including those which exceed 100 mm/km, prior to 6 a.m. local time.
The results on Figure 10 were obtained by combining nine days of data (the days in Table 1 and Table 2, except ID 6) and computing vig (the curve with circles) and vig (the curve with triangles) separately for each two-hour time bin. The periods from 4 to 6 a.m. local time and 6 to 8 p.m. local time (just beyond daytime hours) have much higher vig overbounds (the curve with asterisks) than the periods between 6 to 8 a.m. local time and 4 to 6 p.m. local time. For this reason, we recommend that the definition of daytime for GBAS operations should correspond to the period between 6 a.m. and 6 p.m. local time. Because this is a general guideline based upon data analyzed from throughout Brazil, and because the beginning and end of daytime for ionospheric purposes is most closely related to sunrise and sunset times (which vary throughout the year), more specific guidance is required for the GBAS at GIG. This is addressed further in Section 5. The results in both Table 2 and Figure 11 obtained using the entire Brazil database within 6 a.m. and 6 p.m. local time suggest that the bounding ionospheric contribution to the nominal broadcast vig should be 13 mm/km. Table 2 for Dataset 2 identifies February 25, 2014, as the day with the highest vig when the vig overbound is computed separately for each day, and that highest value is 13 mm/km. In Figure 11, which is a zoomed-in version of Figure 10, when the vig overbound is computed over nine days (including February 25, 2014) but is broken down into two-hour intervals, the maximum vig value within daytime is also 13 mm/km.
Selecting 13 mm/km as the ionospheric component of the nominal vig overbound and performing an RSS with the same 5 mm/km anomalous tropospheric contribution as in CONUS results in a minimum broadcast value of 13.93 mm/km, which would be rounded up to 14 mm/km. It is not known if anomalous tropospheric conditions in low latitudes (e.g., in tropical regions) can be worse than those in CONUS. Large tropospheric delay gradients are caused by weather conditions such as strong cold fronts that create large differences in atmospheric temperature, pressure, and humidity over short distances (van Graas & Zhu, 2011), and these conditions should not be worse in tropical regions (except during hurricanes or cyclones, when aircraft operations would be suspended). Since low-latitude ionospheric spatial decorrelation dominates the RSS of the ionospheric and tropospheric variances, it is unlikely that extreme low-latitude tropospheric conditions are worse by enough to significantly affect this result.

DEFINITION OF LOCAL DAYTIME
As explained above, the definition of "daytime" hours as between 6 a.m. and 6 p.m. local time is really an approximation of local sunrise and sunset times, as ionospheric behavior violating the vig value proposed above appears to be driven by sunrise/sunset rather than specific local times. Therefore, the GBAS installation at Rio de Janeiro should use an approximation of local sunrise and sunset times, which vary throughout the year, instead of specifically using 6 a.m. and 6 p.m. local time. Since the GBAS naturally observes GPS time instead of local time, the daytime transition times need to be referenced back into GPS time to be implemented within the GBAS. Various online almanacs exist for computing sunrise and sunset times for particular locations and times (see "Sunrise and Sunset Times in Rio de Janeiro," (n.d.) and "NOAA Solar Calculator," (n.d.)). Table 3 shows the times of sunrise and sunset ("Daylight") at Rio de Janeiro for the same day (the 21 st ) of each month during the year 2018 based on the detailed tables provided by "Sunrise and Sunset Times in Rio de Janeiro" (n.d.). The trend of shortening and lengthening days by season is evident here. The maximum difference in sunrise times between solstices is shown to be about 1.47 hours (5:04 a.m. "standard" time on December 21 vs. 6:32 a.m. "standard" time on June 21). Note that "standard" time here refers to local time without the "Summer Time" adjustment previously applied in Brazil (prior to 2020) during its summer months. A review of the provided tables in "Sunrise and Sunset Times in Rio de Janeiro" (n.d.), generated for the next 20 years, shows that these times change very little; thus, a daytime conversion table built based on 2018 should be valid for the lifetime of an SLS-4000 at Rio.
This variation in sunrise and sunset times needs to be provided to the GBAS in a convenient format, along with a means for the GBAS to convert between local time, UTC, and GPS time. For now, it will be assumed that the latter is already present within the GBAS. (This will be revisited if a means for conversion needs to be provided). A table of local sunrise and sunset times can be created and inputted into the GBAS at varying levels of precision (e.g., one entry per month), ignoring any shifts caused by local Summer or "Daylight Savings" Time (as these do not affect GPS time). An alternative is to fit polynomial curves to the actual sunrise and sunset times (as a function of day of year) to allow the GBAS to compute estimates of these times for each day. Precise accuracy in computing sunrise and sunset times is not required: Results within ± 5 minutes of the actual times should be sufficient. To the degree possible, these models should be developed offline to minimize the additional calculations required of the GBAS. TA B L E 3 Sunrise and sunset times (in standard time) and length of day at Rio de Janeiro of each month for the same day of the month during 2018 One other consideration is that the days used to establish the minimum nominal vig and the daytime period in this study were all between September and April, corresponding to the more active months of the year for the ionosphere in South America. These months also correspond to days with longer daytime periods, as shown in Table 3. By using the results from these months during the quieter mid-year period (Brazilian winter) when the days are shorter, we may be overly limiting the achievable GBAS performance, meaning that the minimum broadcast vig could be lower and/or that the definition of "daytime" could be extended while still being protected by the minimum broadcast vig and CONUS threat model. Additional data analysis of mid-year days would be needed to verify this hypothesis.

6
TEMPORAL GRADIENT EFFECTS

Temporal gradient estimation using "geometric similarity" method
This section evaluates the temporal gradient effect added to the vig bound on spatial gradients introduced by the time-step method. A new method called "geometric similarity" is proposed to estimate temporal gradients. This method describes differential ionospheric delay computed by the time-step method as a linear combination of spatial and temporal delay differences. If it is possible to extract the component of spatial delay difference from the total delay difference, ionospheric temporal decorrelation can be estimated. The underlying assumption to make this possible is that the ionospheric decorrelation between two points can be measured by both the station-pair method and the time-step method. Figure 12 shows the concept of the "geometric similarity" method with the combination of the time-step method in Figure 12(b) on the right and the station-pair method in Figure 12(a) on the left. 1 1 denotes the ionospheric delay measurement at "Point 1″ within the (assumed) ionospheric thin shell at 1 , and 2 1 is the ionospheric delay measurement at a separate "Point 2″ within the ionospheric thin shell at 1 . 2 2 is the ionospheric delay measurement from "Point 2″ at 2 (the time after Δ from 1 ).
The time-step method was used for estimating ionospheric spatial gradients in Section 4. If we use the symbols described above, the differential delay obtained using the time-step method can be written as (from Figure 12 Assuming that 2 2 and 2 1 are delay estimates at the same point ("Point 2") within the ionospheric thin shell but are at different times with a time interval of Δ = 2 − 1 , the ionospheric delay difference caused by temporal variation can be described as temporal delay dif f erence = 2 2 − 2 1 (4) The differential delay from the time-step method in Equation (3) can be divided into two parts; temporal delay F I G U R E 1 2 Illustration of "geometric similarity" method [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] F I G U R E 1 3 Temporal delay difference estimation procedure [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] difference and spatial delay difference.
Therefore, if we can obtain the spatial ionospheric delay difference ( 2 1 − 1 1 ) using the station-pair method, we can obtain the temporal delay difference ( 2 2 − 2 1 ) as follows: The procedure for calculating the temporal delay differences is shown in Figure 13. First, we search for station pairs which have baseline distances below 50 km for the GBAS application. Then, we select station pairs whose IPP separation vector is almost the same as that traced out by the time-step method (for the same satellite) based on the following criteria. First, the magnitudes of the IPP separation vectors for the station-pair method and the timestep method are expressed as 1 and 2 , respectively. The angle between the two IPP separation vectors from the station-pair method and time-step method is defined as . All station pairs are selected where the absolute difference between magnitudes (| 2 − 1 |) is smaller than 10 meters and the magnitude of the angle ( ) is smaller than 10 degrees. The station pairs that meet the criteria are assumed to be similar enough to support the use of Equation (6) and extract the temporal delay difference. Finally, F I G U R E 1 4 Locations of the GEONET reference receivers in September 2010 [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] the temporal delay difference is divided by Δ to calculate the temporal gradient.

Data used to estimate temporal gradients
The sparsity of the Brazilian GNSS observation networks makes it difficult to use the proposed "geometric similarity" method to quantify temporal effects in the time-step method, since it is difficult to find station pairs which meet the selection criteria described above. Thus, we utilized GNSS Earth Observation Network System (GEONET) data from Japan. GEONET comprises more than 1,200 reference receivers and has a dense GNSS observation network with an average receiver separation of about 20 km (Saito, Suzuki, Yamamoto, Saito, & Chen, 2017) as shown in Figure 14.
Plasma bubbles in equatorial regions often cause large ionospheric gradients, and these events are usually associated with ionospheric scintillation. In the Japanese region, plasma bubbles occur frequently during the equinox season (Nishioka, Saito, & Tsugawa, 2008;Saito, Fujita, & Yoshihara, 2012). Therefore, to obtain statistical information of ionospheric temporal gradients for daytime in Brazil under ionospheric scintillation conditions, we chose daytime data (6 a.m. to 6 p.m. local time) from GEONET during 21 days from September 2-22, 2010. To reflect Brazil's low-latitude ionospheric conditions, we selected station pairs located lower than 25 • N in geomagnetic lat-itude (which is about 33 • N in geographic latitude at the longitudes of Japan).

Results of temporal gradient estimation
Using the database just mentioned, Figure 15(a) on the left shows the temporal delay difference (in meters) as a function of IPP separation distance. Different colors for each temporal delay difference point indicate different sizes of the time window (Δ ) used for the time-step method. Figure 15(b) on the right shows the temporal gradient (in meters per minute), which was computed by dividing the temporal delay difference by the corresponding time window to estimate ionospheric temporal decorrelation over a time interval of one minute as follows: The temporal delay difference increases as the IPP separation distance increases because its corresponding time interval between the ionosphere states at 1 and at 2 gets longer. When the temporal delay difference was divided by the time window over which it was estimated, the temporal gradients appear to be relatively consistent regardless of the IPP separation distance.
As explained above, we estimated temporal gradients over a time interval of one minute using GEONET data on 21 days in September 2010. The sample standard deviation of temporal gradient estimates is 7.218 mm/min. When we consider IPP separation distances of about 4-6 km (a typical range corresponding to a one-minute time window), the standard deviation of 7.218 mm/min can be converted to 1.2∼1.8 mm/km = 7.218 mm/4∼6 km. Therefore, the resulting temporal effect on estimated spatial gradients is about (or somewhat below) 2 mm/km when a time window of one minute is used for the time-step method.
Because this is only an approximation of the temporal gradient effect, it is not practical to remove this effect from the overbounding vig derived in Section 4. What these estimates of temporal gradients using the proposed method indicate is that vig estimates derived using the time-step method that necessarily include a temporal gradient component are somewhat conservative but not too conservative for the Brazil application, since the estimated temporal component of ∼ 2 mm/km is a small fraction of the suggested bounding vig of 13 mm/km.

SUMMARY AND REMAINING STEPS
This paper has demonstrated that, while the value of vig that bounds nominal ionospheric spatial gradients has been determined to be approximately 4 mm/km in midlatitude CAT I GBAS operations, a significantly higher value should be broadcast for vig by the GBAS facility at Rio de Janeiro (GIG) to bound low-latitude ionospheric conditions. The results in Section 4 demonstrate that a minimum broadcast vig value of 14 mm/km (including an ionospheric component of 13 mm/km and a smaller tropospheric component) is sufficient to bound nominal ionospheric behavior during daytime. "Daytime" conditions have been shown to generally apply to the period between 6 a.m. and 6 p.m. local time, or, more particularly, between local sunrise and local sunset. This paper also evaluates the effect of ionospheric temporal gradients included in the estimates of ionospheric spatial gradients computed using the time-step method. With the proposed "geometric similarity" method, the standard deviation of temporal gradients over a oneminute interval is estimated to be about 2 mm/km, which is a small fraction of the overall bounding value of 13 mm/km. Because the times of local sunrise and sunset vary both in terms of local time and GPS time (which is observed by the GBAS), a practical means is needed to allow the GBAS at Rio de Janeiro to determine when to switch to and from "daytime" conditions (and thus to activate or deactivate GBAS approach service). An approach to do this is proposed in Section 5 but needs review and further consideration. A key issue is how often these switching times (in terms of GPS time) need to be updated over the course of one year to compensate for variations in local sunrise and sunset.
These results should be sufficient to support the approval and certification of CAT I precision approach operations using the existing GBAS at Rio de Janeiro during the daytime hours defined here. In future work, these analysis tools can be used to determine how to extend GBAS operations into the early morning and early evening periods while maintaining integrity against the moresevere ionospheric conditions that exist outside daytime hours. A re-examination of the seasonality of the gradients that motivate the proposed vig bound is recommended, as it may be possible to lower this bound and/or extend the definition of daytime hours during the middle of the year, when large gradients are generally not observed.

A C K N O W L E D G E M E N T S
This paper is a revised version of the paper previously published in the 2019 International Technical Meeting of The Institute of Navigation proceedings.