Leveraging Geodetic GPS Receivers for Ionospheric Scintillation Science

We demonstrate scintillation analysis from a network of geodetic Global Positioning System (GPS) receivers which provide data at 1‐second resolution. We introduce proxy phase (σTEC) and amplitude (SNR4) scintillation indices and validate them against the rate of change of TEC index (ROTI) and S4. Additionally, we validate scintillation observations against a Connected Autonomous Space Environment Sensor scintillation receiver. We develop receiver‐dependent scintillation event thresholding using hardware‐dependent noise variance. We analyze 6 days adjacent to the 7–8 September 2017 geomagnetic storm, using 169 receivers covering magnetic latitudes between 15° and 65° in the American longitude sector. We leverage the available spatial sampling coverage to construct 2‐D maps of scintillation and present episodic evolution of scintillation intensifications during the storm. We show that low‐latitude and high‐latitude scintillation morphology match well‐established scintillation climatology patterns. At midlatitudes, spatiotemporal evolution of scintillation partially agrees with known scintillation patterns. Additionally, the results reveal previously undocumented midlatitude scintillation‐producing structures. The results provide an unprecedented view into the spatiotemporal development of scintillation‐producing plasma irregularities and provide a resource to further exploit scintillation evolution at large spatial scales.


Introduction
The ionosphere significantly alters traversing radio signals at frequencies <3 GHz (P.531-14, 2019), which include Global Navigation Satellite Systems (GNSS) signals. In particular, small-scale (<1 km) plasma density irregularities affect the L-band radio signals by means of Fresnel diffraction, causing the signal's amplitude and phase to rapidly fluctuate in a stochastic manner. A theoretical description of the diffracting scatter has been reviewed in detail (Kintner et al., 2007;Priyadarshi, 2015;Yeh & Liu, 1982). On the other hand, larger-scale plasma density irregularities (>1 km) impose fluctuations in received phase only, by means of signal refraction.
Extensive surveys of occurrence and climatology of both phase and amplitude scintillation have been conducted (Aarons, 1982;Alfonsi et al., 2011;Basu et al., 1988;Jiao & Morton, 2015;Wernik et al., 2003), primarily focusing on the low-and high-latitude regions. It has often been considered that necessary conditions for ionospheric density irregularities are confined to the equatorial and low-latitude regions (by Rayleigh-Taylor instability) and high latitudes (convection-driven instabilities and particle precipitation). The midlatitude ionosphere has not generally been perceived as a possible space weather threat for radio signal scintillation, with exceptions of equatorward expansion of high-latitude convection-related dynamics (Aarons, 1982;Kintner et al., 2007). Even though significant amplitude scintillation was observed from the upstate New York during a geomagnetic storm (Ledvina et al., 2002), such scintillation occurrence geolocation and timing are not predicted by conventional scintillation models. For example, empirical scintillation models, wideband ionospheric scintillation model (WBMOD) (Secan et al., 1995(Secan et al., , 1997, and the global ionospheric propagation model (GIM) (Béniguel, 2002) predict no signal impairments at midlatitudes.
The historical surveys and existing models arise from insufficient availability of ground-based infrastructure for ionospheric scintillation monitoring. Many of the existing scintillation receivers are located at low and high latitudes, and only a few scintillation receivers operate at midlatitudes. Hence, a comprehensive analysis of ionospheric scintillation at midlatitudes, by means of conventional amplitude (S 4 ) and phase ( ), has been impossible. Generally, the scintillation receivers cover localized areas; hence, the climatological statistics and occurrence timing is based on localized observations. We utilize a 1-Hz GNSS receiver network operated by UNAVCO, which covers an area from low (15 • magnetic latitude [MLAT]) to high latitudes (65 • MLAT) over North America with data availability spanning back to 2011. Therefore, the network analysis allows the creation of instantaneous two-dimensional (2-D) maps of scintillation at unprecedented scales. The network consists of hundreds of receivers (∼200 in 2011 and ∼850 in 2020), where the vast majority of the receivers monitor Global Positioning System (GPS) signals only; thus, we limit our study to the GPS.
Utilization of high-rate geodetic receivers is an attractive and potentially ground-breaking avenue for augmentation of GNSS scintillation science. Recently, several studies demonstrated the capability of high-rate geodetic receivers on a single-receiver comparison cases with colocated ionospheric scintillation monitors (Juan et al., 2017;Luo et al., 2020;Nguyen et al., 2019). The preliminary results were focused on the equatorial scintillation. Further, Luo et al. (2020) extended the capability to image amplitude scintillation over China. Contrary to the other studies, we do not attempt to calibrate the values of our scintillation indices to S 4 and , since we have no colocated ionospheric scintillation monitor data available. We build upon the recent results and present scintillation signal processing for a diverse receiver network, capable of imaging both phase and amplitude scintillation and monitoring instantaneous scintillation patterns over a large spatial area. The purpose of this work is to promote the utility of 1-Hz (hereafter referred as high-rate) geodetic receivers for scintillation science; specifically, it is a diagnostic tool that provides an unprecedented insight into small-scale dynamics of large-scale density features, such as storm time midlatitude dynamics (Ledvina et al., 2002;Mrak et al., 2020).
We derive and validate alternative scintillation indices based on the signal-to-noise ratio (SNR) and total electron content (TEC). Although each parameter is affected by different physical phenomena (diffraction and refraction, respectively) (McCaffrey & Jayachandran, 2019), we do not attempt to distinguish among them. We introduce and validate receiver hardware-dependent thresholding for the scintillation event decision. We demonstrate that the introduced indices are linearly correlated with S 4 and ROTI. Comparison with a Connected Autonomous Space Environment Sensor (CASES) scintillation receiver from Dallas, TX, further supports the validity of the introduced indices. The scintillation events are used to produce instantaneous scintillation maps. We demonstrate the utility of the scintillation maps for spatiotemporal analysis of evolving scintillation. We analyze 6 days' worth of data adjacent to 7-8 September 2017 geomagnetic storm. The low-and high-latitude portion of the receiver network observed typical low-and high-latitude scintillation patterns. Additionally, the storm had a striking impact on the midlatitude GPS receivers, exhibiting both amplitude and phase scintillations.

Methodology
We utilize publicly available data from the high-rate UNAVCO GPS receiver network, with available receiver distribution on 8 September 2017 depicted in Figure 1. There is a heavy sampling bias toward the west coast of the continental United States (CONUS). The total distribution of receiver hardware types on that day is presented in Figure 1b. The receiver network covers area between 15 • and 65 • MLAT in the northern hemisphere. Most of the receivers are placed within the contiguous U.S. longitude span, whereas the high-latitude coverage comes from the Alaskan sector. In this study, we define midlatitudes as the area between 30 • ≤ MLAT ≤ 60 • in the Northern Hemisphere.
Normally, scintillation studies use specialized GNSS scintillation receivers that provide amplitude (or power) and carrier phase at subsecond (usually 50 Hz) data rate. These measurements are used to calculate Φ and S 4 scintillation indices. The high-rate geodetic receivers do not provide signal amplitude, and the carrier phase is contaminated by receiver-imposed interference. We chose to utilize SNR provided by the Receiver Independent Exchange (RINEX) files, as a substitute for signal amplitude. SNR has an inherently bigger variance than the amplitude, as it is influenced by the variance of broadband noise. Additionally, SNR information provided by the RINEX files is by definition (computed) receiver dependent, and it is defined as an SNR of the demodulated carrier, namely, the carrier-to-noise ratio. We build upon successful demonstrations of SNR utilization for scintillation studies (Rodrigues & Moraes, 2019;Thompson et al., 2008), and we evaluate the behavior of scintillation indices for different receiver types.
We found an impediment with Trimble NETRS carrier phase measurements, which is depicted in Figure 2. This example represents the general trend bound to this particular receiver make, not to a specific receiver This panel depicts numerous abrupt jumps in the carrier phase. The consequence of this impairment is that we can't use conventional phase scintillation index Φ as it is defined on L1 carrier phase measurements. The peculiar nature of the carrier phase and its use for scintillation studies has been discussed by Beach (2006). A carrier phase combination of L2P and L1CA channels in the bottom two panels shows that the abrupt changes cancel out. The resulting detrended phase combination Δ(L2P-L1CA) in panel (c) and its high-pass filtered derivative show receiver hardware-independent carrier phase perturbations, which can be used for extracting ionosphere scintillation effects. Trimble NETRS receivers comprise more than half of the UNAVCO receiver network, and thus, we adopt the use of carrier phase combination instead of the L1 carrier phase throughout the network. As we show in the next section, carrier phase combination is proportional to the TEC, and thus, we utilize this property and the TEC as a measure of phase fluctuations. As suggested by Beach (2006), the use of a TEC, that is, a weighted carrier phase combination, can eliminate the receiver-imposed errors to phase measurements. Measurements sampled at 1 Hz tend to undersample scintillations as the scintillation spectra extend to subsecond time scales. The carrier phase is not an issue as the power law dependence shall continuously extend into a larger spatial scale (slower time scale at the receiver). It has been demonstrated that a 1-Hz receiver and a colocated scintillation receiver measure morphologically the same scintillation spectra above 0.5 Hz (Béniguel et al., 2009). On the other hand, amplitude scintillation spectra has a low-frequency cut-off (phase screen approximation) at the Fresnel scale r F (r F = √ 2 Z, Z being the distance between receiver and irregularity) (Kintner et al., 2007), which translates into Fresnel frequency f F = r F v 0 , where v 0 is the irregularity drift. Thus, a 1-Hz receiver with a Nyquist frequency at 0.5 Hz oversamples only scintillation-producing irregularities with effective drift velocity v 0 ≤ 180 m/s (Rino, 1979). This calculation assumes isotropic irregularities at Z = 350 km, with v 0 perpendicular to the line of sight.

Data Processing
We utilize GPS data from RINEX files where we take SNR at L1 ( 1 = 1575.42 MHz) and deduct slant TEC (sTEC) from carrier phases L1CA and L2P ( 2 = 1227.6 MHz) expressed as ranges where cycle slips are accounted for with a method of Blewitt (1990), carrier phase ambiguities with a method from Coster et al. (1992), and 1 TECu = 10 16 electrons per meter squared. The sTEC is derived from the phase accumulation property that the phase advance is inversely proportional to f 2 and proportional to sTEC ( L= 40.3 2 sTEC) (Mrak et al., 2018). The sTEC is then converted to vertical (hereafter referred as TEC) assuming the thin shell ionosphere approximation via mapping function F(Θ) (Klobuchar, 1987) (2) via zero TEC method (Rideout & Coster, 2006). Here, we assume the height of the ionospheric piercing point (h IPP ) to be 350 km, R e is the radius of the Earth, and Θ is the elevation angle. We apply a (sixth order) high-pass filter, with 0.1-Hz cut-off frequency and Butterworth response (Fremouw et al., 1978) to the sTEC and SNR, obtaining TEC and SNR, respectively. Scintillation indices TEC and SNR 4 are then computed with a running 60-second standard deviation filter, normally used to compute phase scintillation index .
We further adjust low elevation measurements due to oblique angle propagation through the irregularity layer using the approach by Alfonsi et al. (2011) and Spogli et al. (2009). Then, SNR 4 can be expressed as Due to the problem with the carrier phase mentioned above, we compare the phase scintillation index TEC against the rate of change of TEC (ROT) index (ROTI) similar to the early work of Beach and Kintner (1999). We compute ROTI, utilizing a running 60-second standard deviation filter, as defined by Pi et al. (1997): where the ROT is differential TEC computed at 1-second resolution. Furthermore, we estimate equivalent S 4 by utilizing the SNR measurements, building upon promising results of recent case studies (Luo et al., 2020;Rodrigues & Moraes, 2019). We convert SNR into linear units of intensity I = 10 SNR∕10 and then compute the index on a running 60-second window. Finally, we account for the elevation angle as introduced in Equation 5.  Kintner et al. (2007) discussed signal intensity-SNR relationship, where SNR is averaged over a period of 1 second. The exact formulation used by receiver types used in this study is unknown. On one hand, the introduced technique for S 4 estimation has been proven to correlate well with conventional S 4 computed directly from the signal intensity (Luo et al., 2020;Rodrigues & Moraes, 2019;Thompson et al., 2008). On the other hand, the primary drawback of using SNR is a fade smearing due to the temporal averaging, as demonstrated by Jiao et al. (2016).

Event Definition
We have discussed data impairments and mitigation of geodetic receivers used in the UNAVCO network.
Here we discuss the "scintillation event" selection procedure, which takes into account inherent noise level variability due to diverse receiver hardware selection. In the event selection criteria, we build upon the statistical survey of Jiao and Morton (2015), with additional receiver-specific factor. Figure 3 presents a distribution of receiver hardware setups and related histogram of noise levels of the introduced indices. We chose to use a daily median value of phaseT EC and amplitudeŜNR 4 scintillation index as a measure of the receiver noise floor. We compute the statistical distribution of these levels in a full month of January 2018, excluding days with planetary K (Kp) index Kp ≥ 4 (14 and 24 January). All receivers were located in the Northern Hemisphere, where the total background TEC is the lowest on a yearly basis. Figure 3a depicts average hardware distribution in the time period of this analysis. Figure 3b shows a distribution of receiver noise levels among all available receivers. Figure 3c breaks down the receiver noise levels for each receiver setup. Overall, bothT EC andŜNR 4 have a large spread with a factor of ∼5 in TEC and a factor of ∼10 in SNR 4 . While the degree of spread in TEC is similar among different receiver setups, with the worst performer being Trimble NETR8, Septentrio PolaRX5 is notably a receiver with the smallest noise level and variance in SNR.
To cope with receiver-dependent noise levels, we define a threshold parameter T r , for ∈ [ TEC , SNR 4 ], for each receiver r computed on a daily basis. The threshold values are defined as where a value of 2.5 is a fixed standoff distance from the receiver noise level̂r. Although this number is empirically chosen, as demonstrated below, this threshold effectively separates real scintillation signal and noise. Scintillation events are then defined separately for TEC and SNR 4 as follows: running median (60-second length) value of a scintillation index has to continuously exceed the computed threshold (Equation 8) for a minimum duration of 2 minutes. Additionally, multiple events with temporal separation shorter than 5 minutes are merged together. The introduction of the variable threshold value T r depending on the noise level is a modification from a fixed cut-off method used by Jiao and Morton (2015), because a fixed threshold creates a bias between different types of receivers.
The event selection is pictorially demonstrated in Figure  . The bottom panel in Figure 4 demonstrates the event selection result by the thick line. The increased SNR 4 due to the system glitch was rejected by the minimum length criterion. On the contrary, the event selection procedure flagged elevated SNR 4 between 00:45 and 1:00 UT as a scintillation event.
The second example in Figure 4b depicts a case of a longer phase scintillation event (labeled as "1") colocated with a TEC gradient at 0:40-1:30 UT. Additionally, there was a secondary, more subtle enhancement in the TEC around 2:00 UT (marked as "2"). This scintillation event consists of a period where TEC was below the threshold (in the middle of the event). This part was included by the minimum separation criterion, which merges individual scintillation events with a separation shorter than 5 minutes into a single event.

Case Study
We analyze and validate the introduced indices and scintillation event selection on an event study. We have processed 6 days' worth of data surrounding the 7-8 September 2017 geomagnetic storm. The solar wind and geomagnetic indices for this period (5 to 11 September 2017) are presented in Figure 5. The solar wind data show that two consecutive shocks hit the magnetopause on 7 September, with a time separation of about 24 hours. While the first shock arrived with a predominantly northward interplanetary magnetic field (IMF), the latter one was accompanied by a strong (∼ − 30 nT) southward IMF. The latter shock facilitated the enhancement of the ring current (SYM/H index) and increased high-latitude geomagnetic activity measured by the auroral electrojet (AE) index. There were two episodic AE intensifications on 8 September; both were related to storm development (negative excursion of the SYM/H). Note. Locations are given in geographic longitude (GLON) and geographic latitude (GLAT).
The storm was particularly intriguing, as it severely perturbed the ionosphere at the longitude sector of northern America, and thus, it was well sampled by the UNAVCO GPS network. The ionosphere exhibited several distinct perturbations at this local time sector, in the form of equatorial plasma bubbles (EPBs), severe auroral activity, and multiple TEC gradients at midlatitudes (over CONUS). Impulsive perturbations have been reported (Aa et al., 2019;Mrak et al., 2020;Zakharenkova & Cherniak, 2020), but small-scale density irregularity distribution remains unknown, due to the lack of scintillation receivers. Severe space weather impacts on GNSS over the contiguous United States have been presented (Yang et al., 2020); however, associated scintillation or small-scale irregularities have not been yet reported. We take advantage of the UNAVCO GPS receivers and demonstrate the scintillation processing for this event.

Evaluation of Scintillation Indices and Event Selection
We analyze scintillation indices computed from a subset of six receivers during the time frame of the 7-8 September storm main phase. The chosen receivers are listed in Table 1, and we focus only on the three most affected lines of sight for each of the receivers. The corresponding trajectories of ionospheric pierce points (IPPs) at 350-km altitude are depicted in Figure 6 for the period shown in Figure 7. Derived TEC and scintillation indices for these receivers are presented in Figure 7. Receivers are chosen in a way to cover a large span in longitude and latitude and to present time series plots of the three most frequent hardware setups. Figure 6. Locations of receivers from Table 1, with trajectories of ionospheric piercing points from the chosen satellites.

10.1029/2020RS007131
Figure 7. Example measurements from three receivers, introduced in Figure 6. Each panel consist of (top to bottom) estimated vertical TEC (vTEC), ROTI, σ TEC , SNR 4 , and S 4 . Each receiver represents different hardware setup, presented in Table 1. Red dashed lines are the associated thresholds T r .

10.1029/2020RS007131
Each panel in Figure 7 consists of computed TEC (Equation 1), ROTI, and SNR-derived S 4 , TEC , and SNR 4 . All presented observations were masked with a 30 • elevation cut-off. The TEC plots show dramatic TEC perturbations, some exceeding 10 TECu at time scales of a few minutes. In general, ROTI values correlate well with these gradients, revealing the existence of irregularities with temporal scales in the order of 1 second. A comparison between ROTI and TEC indicates a linear correlation, showing no morphological differences regardless of the receiver type and intensity of irregularities. While most of the receivers observed elevated TEC perturbations by virtue of both indices, the event threshold T TEC ≈ 0.01 (translates toT EC ≈ 4 ⋅ 10 −3 ) for all receivers, regardless of the receiver type.
The bottom two panels of each receiver in Figure 7 serve as a comparison between amplitude scintillation indices. Like in the TEC case, the indices are visually well correlated but show strikingly different variance (noise level) among different receivers. A large span in the thresholds T SNR 4 is expected based on the preliminary quiet day analysis (Figure 3). In the given example, Septentrio PolaRx5 receivers (wmok and hdil) have the threshold level of a factor of ∼2 smaller than Trimble receivers. We find that the introduced SNR 4 index has a bigger dynamic range, due to the fact that it is not normalized.
Another striking observation is enhanced amplitude scintillation, as four of the receivers were located within CONUS: Washington state (p413), Oklahoma (wmok), Illinois (hdil), and California (p209). Less strikingly, receivers at lower latitudes-Jamaica (cn12) and Mexico (oxum)-were as well subject to enhanced scintillation. Specifically, receiver oxum recorded extreme amplitude scintillation with S 4 ≥ 0.5, accompanied by numerous loss of locks, which prevented the calculation of TEC index. Phase scintillation was colocated with amplitude scintillation at all instances. Receiver p413 measured only increased phase scintillation, indicating that it got affected by different kinds of ionospheric structures. Lastly, enhanced amplitude scintillation is well correlated with phase scintillation TEC and ROTI, generally colocated with the steep TEC gradients.
The correlation between ROTI and TEC observed with all receivers is not surprising. As we introduced in the previous section via Equation 1, they should be linearly correlated. The advantages of TEC over ROTI are a deterministic amplitude and phase response by the prescribed filtering operation and a straightforward connection with the TEC power spectral density (PSD TEC ) following the morphology derived for phase scintillation (Béniguel et al., 2009;Secan et al., 1995), where the low limit in the integration is a chosen frequency cut-off at 0.1 Hz. Correlation between SNR 4 and TEC does exists from cases of measured amplitude scintillation. The latter observation is also expected as several studies demonstrated correlation and casual relationship between ROTI and S 4 (Carrano et al., 2019;Liu & Radicella, 2019;Yang & Liu, 2016). Correlation deviation in this formalism is thought to be due to viewing geometry and irregularity drift velocity (Carrano et al., 2019;Liu & Radicella, 2019).
The importance of the hardware-dependent scintillation event classification is bolstered by means of large-scale statistical relationships: TEC -ROTI and SNR 4 -S 4 . A mutual comparison between the indices is presented as 2-D histograms in Figure 8. The histograms encompass data from the receivers in Table 1 on a period of 48 hours (7 and 8 September 2017), with a total of ∼4.5 million data points. The histograms show two major populations: a correlated group clustered along a linear correlation line and the second uncorrelated group along both axes. The latter group exists due to receiver impairments such as cycle slips and processing artifacts.
We then apply the scintillation event decision criteria and redo the correlation analysis in Figure 9. Panel (a) shows the results from the same data set as in Figure 8. The criteria effectively rejected the uncorrelated groups, and the results show a linear relationship between both TEC -ROTI and SNR 4 -S 4 . Linear correlation with a high correlation coefficient indicates that the introduced scintillation indices are adequate substitutes for ROTI and S 4 . The bottom panel (b) presents analysis taken from all receivers in Figure 11 processed in this case study. The linear correlation line is plotted on each panel, including the line parameters and the correlation coefficient R. We find the correlation coefficient between TEC and ROTI R = 0.9 and SNR 4 and S 4 R = 0.83.

Comparison Against CASES Scintillation Receiver
We compare measurements from the wmok receiver located in Oklahoma with a CASES scintillation receiver (Crowley et al., 2011) deployed at UT Dallas, TX. The receivers were spaced ∼240 km apart. As the UT Dallas receiver was operating only until 2:30 UT on 8 September, we utilize data recorded in this period. The CASES receivers sample carrier phase and signal intensity of the L1CA channel at 100-Hz resolution and output scintillation indices and S 4 every 100 seconds. Thereby, the rate deviates from the conventional approach (60 seconds), which we adopted in this work. For a purpose of this comparison only, we used the same 100-second standard deviation window for the wmok receiver. CASES receivers provide estimates of the TEC and SNR at 1-second rate.  We analyze observations from three satellite links, where we computed the introduced scintillation indices ( TEC , SNR 4 ) at a 100-second rate for wmok receiver, which we compare with the CASES-provided and S 4 . The comparison is presented in Figure 10, where the top row consists of frames with observations taken by the CASES receiver and the bottom row of the wmok. All observations were masked with 30 • elevation angle cut-off; the elevation angle is depicted alongside the wmok-TEC plots. Estimated TEC is plotted in the top panel, whereas the derived TEC is in the second panel. It is plotted side by side with the CASES-provided phase scintillation index . All three lines of sight show that TEC follows the trend and dynamics of . The bottom two panels show the SNR in the third panel, whereas the SNR 4 and S 4 indices are in the last panel. The latter amplitude scintillation indices have large variance owing to the receiver's performance, as indicated by large variations in the SNR panels. The SNR 4 index was derived from the SNR and follows the trend of the S 4 which is provided as a standard CASES output.
The general trend in the scintillation indices (red lines) is morphologically similar between both receivers. While the middle frames, monitoring satellite G24, show almost identical trends in both TEC and SNR 4 ; the other links have a distinct difference. The difference originates in different electron density structure as measured by the TEC. Notable differences can be observed between satellites G10 and G32, where both receiver measured double-slope gradients but at the opposite links. The difference can be ascribed to the 240-km distance between the receivers. Nevertheless, both receivers measured increased amplitude scintillation in the plateau region between the density gradients. This was measured by the CASES receiver in the G10 link at ∼1:15 UT and by the wmok at about the same time. It should be noted that despite that increased variance in raw SNR data, namely, amplitude scintillation, is clearly observed, the value of the scintillation index S 4 is below 0.2 usually used as a scintillation detection threshold (Béniguel et al., 2009;Jiao & Morton, 2015).
In aggregate, the CASES scintillation receiver provides independent validation for the introduced phase scintillation index substitute, TEC , that follows the trend of . Additionally, CASES-derived amplitude scintillation index S 4 measured instances of elevated amplitude scintillation, which were also indicated by the SNR-derived SNR 4 . In comparison with the closest UNAVCO receiver, it is found that both receivers measured the same scintillating structures, both with phase TEC and amplitude SNR 4 scintillation indices. Figure 11. Same format as Figure 1, with spatial distribution of receivers chosen (total of 169) to produce irregularity maps in this case study.

Scintillation Maps and Event Analysis
The spatiotemporal evolution of the scintillation indices during the storm is obtained by virtue of scintillation maps from chosen receivers presented in Figure 11. A total number of 169 receivers were selected based on the criteria to reduce receiver density to one receiver per 2 • GLON × 2 • GLAT bin. The total number and specific receiver hardware contribution to the subset of the receiver network is presented in Figure 11b. We present the irregularity maps in Figure 12, covering the 7-8 September storm at a time resolution of 1 hour. Irregularity maps consist of TEC and SNR 4 indices, which are compared to the ROTI maps at each epoch. Each frame consists of data points collected in the first 5 minutes after the image time stamp.
The maps in Figure 12 show how well the regions of elevated ROTI correlate with scintillation occurrence. Each epoch consists of two maps, where ROTI maps are the top panel and the scintillation maps are the bottom one. All measurements are mapped to 350-km latitude and superposed on top TEC maps. While we plot all the ROTI observations, the scintillation indices underwent the classification process described in section 2.2. At the beginning of the storm's main phase in Figure 12a, the area with high ROTI correlated with phase scintillation TEC . The disturbance lied within a storm-enhanced density plume (area of enhanced TEC), extending from the central United States up to northern Alaska. In the next frame (panel b), the morphology of scintillation at the poleward portion of the CONUS remained unchanged, while a meridional density depletion (seen in the TEC map) over the eastern United States hosted another disturbance with elevated ROTI and both amplitude and phase scintillation. For the evolution of the TEC depletion, we refer readers to recent storm studies Zakharenkova & Cherniak, 2020).
An hour later in Figure 12c, the maps show a change in scintillation morphology in the vicinity of the midlatitude trough (marked in panel a) located near the United States-Canada border. There were two distinct meridionally separated scintillation regions: one poleward in the midlatitude trough (within the auroral oval) and the other one equatorward of this trough. The scintillation elongated along the secondary density trough persisted through 4:30 UT. Additionally, there were numerous instances of observed amplitude scintillation within the secondary density trough. These observations were previously described in time series plots of receiver p413 in Figure 7d. The scintillation lingered over the CONUS area for at least another hour (panel d), before they began to dissipate away.
Another region of intense amplitude scintillation emerged at about 1:30 UT (Figure 12c), located within the high-density area near 20 • MLAT, just west of 0 • magnetic longitude meridian. Amplitude scintillation, accompanied with phase scintillation (in conjunction with elevated ROTI), persisted until the end of the event analysis at 4:30 UT. The scintillation observations were located deep within the high-density (the equatorial ionization anomaly) area. Considering the timing aspect of it (premidnight local time) strongly suggests that the scintillation was a consequence of EPBs (e.g., Béniguel et al., 2009). At the same time, high-latitude activity was strictly confined within the auroral oval, over the Alaskan sector.
Lastly, we analyze total scintillation occurrence during 6 days adjacent to the storm under investigation. We present scintillation occurrence as a function of time in Figure 13. We took the same data as used to create the instantaneous maps of scintillation discussed above and reduced the spatial dimensions into normalized scintillation parameters. The ROTI time series plot in panel (a) represents the median (< ⋅ >) value of ROTI at a given epoch. On the contrary, the sole use of this statistical moment would not comprehensively represent the scintillation occurrence due to the fact that the total number of scintillation events varies with time as defined in section 2.2. Therefore, we multiply the median value of a scintillation index with the total number of events N recorded at a given epoch. Thus, the total parameter encompasses information about (median) scintillation intensity and spatial spread of it. Resulting normalized phase scintillation < TEC > ⋅N and amplitude scintillation <SNR 4 > ⋅N are plotted in the bottom two panels (respectively).
The phase scintillation occurrence in Figure 13b has three distinct peaks (enumerated I-III), which correlate well with auroral electrojet (AE index) intensifications in the background (thin blue line). The correlation is expected as the AE is a proxy measure of auroral activity which is well established as the predominant source of phase scintillation at high latitudes. The second (II) increase in TEC , however, is disproportional to the other two instances (I and II). This increase occurred in conjunction with the extreme increase in amplitude scintillation SNR 4 , as shown in Figure 13c. This intensification occurred at the dip of the geomagnetic storm, indicated by the negative deflection in SYM/H index (shown in the bottom panel). Finally, this storm time increase correlates well with elevated ROTI; however, the other two instances of elevated phase scintillation index TEC do not show a contemporary increase in the ROTI. This is due to the fact that only a localized area of the receiver network coverage experienced the scintillation; hence, the median operator did not pick up the increases in ROTI, which did not undergo the scintillation event classification decision process.

Summary and Conclusions
We have introduced an alternative method to obtain ionospheric scintillation indices from geodetic GPS receivers with a 1-Hz data rate. We have discussed limitations imposed by opportunistic data source, which has limited temporal sampling range. We have introduced and demonstrated the efficiency of a hardware-dependent scintillation event classification by virtue of direct comparison with ROTI and SNR-derived S 4 . We showed that the introduced amplitude scintillation SNR 4 index is more sensitive to weak events (cf. Figure 7) than the S 4 , due to the fact that SNR 4 is not normalized. The introduced phase scintillation index TEC shows a linear correlation with ROTI, with a correlation coefficient of 0.90. Another property of the introduced scintillation indices is equal signal processing treatment, where the filtering function has a deterministic impulse response, in contrary to the ROT and ROTI.
The introduced processing was applied to the UNAVCO GPS data set; we leveraged its large spatial coverage to produce scintillation maps. The most profound virtue of large spatial coverage is the ability to examine scintillation evolution at a continental scale covering the longitude sector of the CONUS. We demonstrated the potential importance of the data product with a case study of the 7-8 September 2017 storm. The results reveal episodic scintillation occurrence and spatiotemporal evolution in an area covering over 50 • MLAT, considered to be primarily within the midlatitude ionosphere. Six days of scintillation observations adjacent to the 8 September storm show a good correlation between increased auroral activity and the phase scintillation index TEC . A disproportional increase in both phase and amplitude scintillation was observed during the storm main phase. The spatiotemporal evolution of the scintillation geolocation was analyzed by virtue of scintillation maps, namely, high-latitude scintillation was predominantly confined to the area within the auroral oval and in the vicinity of the trough equatorward boundary, characterized by an elevated phase scintillation index. Because this region lacks amplitude scintillation, the elevated phase scintillation index is likely due to phase variations associated with fast-moving density structures. Low-latitude scintillation was predominantly confined within the equatorial ionization anomaly at the premidnight local time sector, consisting of both elevated amplitude and phase scintillation, a morphology of Fresnel scatter. Lastly, additional large-scale density trough was located predominantly at midlatitudes and was associated with elevated amplitude and phase scintillation. The latter finding is a novel observation, made available by the utilization of geodetic receivers.
The presented case of storm time spatial evolution of GPS scintillation at midlatitudes shed a new light on past observations of scintillation from upstate New York (Ledvina et al., 2002). Additionally, the promising results of scintillation event classification and demonstrated utility of scintillation maps allow a comprehensive retrospective analysis of thus far ignored scintillation occurrence at midlatitudes. The data availability goes back to 2011, as a statistical survey is currently underway. While scintillation occurrence at low and high latitudes agrees well with established scintillation climatology (Aarons, 1982;Basu et al., 1988;Béniguel et al., 2009;Kintner et al., 2007;Secan et al., 1995Secan et al., , 1997, the preliminary results from receivers at midlatitudes partially agree with historical morphology of midlatitude scintillation (Aarons, 1982;Kintner et al., 2007). That is, an increase of phase scintillation is predominantly due to TEC gradients in the vicinity of the trough equatorward boundary and increased plasma convection within the trough. The results of the presented case study indicate that there are other midlatitude mechanisms producing both amplitude and phase scintillation. The climatology and controlling parameters of these events are mysterious, as it appears that they occur during large storms (Aa et al., 2019;Ledvina et al., 2002;Zakharenkova & Cherniak, 2020). Leveraging 1-Hz geodetic receivers, such as the one operated by UNAVCO, could be utilized to gain insight into the midlatitude scintillation phenomenon through a comprehensive retrospective analysis.