A Method for Automatic Detection of Plasma Depletions by Using GNSS Measurements

Enhanced scintillation activities observed at transionospheric radio signals are often correlated with slant total electron content (STEC) depletions in the equatorial ionosphere. In this study, the data derived from high‐rate Global Navigation Satellite System (GNSS) receivers were used to analyze the observed STEC depletions, commonly associated with plasma bubbles causing radio scintillations in the equatorial ionosphere. We found that the observed STEC depletion can be described by a wedge‐shaped structure. To quantitatively describe the structure of STEC depletions, we developed an effective method to routinely characterize the depth and the width of equatorial plasma depletion in automatized data screening. The developed method has been validated by analyzing data obtained from mostly African GNSS stations in 2014 and 2015. The results confirm current knowledge regarding the seasonal occurrence of radio scintillations and related bubbles. The detection results are compared with those from other published plasma bubble detection techniques.


Introduction
The equatorial ionosphere with its complex dispersive features is of special interest in space weather research, which supports the operation of space-based radio systems such as the Global Navigation Satellite System (GNSS) and remote sensing radars. The equatorial ionosphere is commonly characterized by the occurrence of enhanced GNSS signal scintillations in the evening hours (e.g., Hlubek et al., 2014;Kriegel et al., 2017). As has already been described many years ago, for example, by DasGupta et al. (1983), the enhancement of amplitude scintillations is inherently associated with total electron content (TEC) depletions. Seemala and Valladares (2011) presented the results of a statistical study by applying an automatic detection algorithm to TEC depletions observed from the Low-Latitude Ionospheric Network, operated in South America during 2008. The method uses a band-pass filter technique applied to calibrated vertical TEC data to detect TEC depletions. To provide statistical information about plasma depletions at low latitudes, Magdaleno et al., (2012Magdaleno et al., ( , 2013 developed an ionospheric plasma bubbles seeker, including an automatized detection algorithm, in Java application. Nishioka et al. (2008) used the difference between the daytime rate of TEC index (ROTI) and the evening ROTI for each day at each GPS station to detect plasma depletions or bubbles during the years 2000. Recently, Blanch et al. (2018 improved the characterization and modeling of plasma depletion by modifying a previously developed technique for detecting medium-scale traveling ionospheric disturbances (Hernández-Pajares et al., 2006). Based on former scintillation studies by Hlubek et al. (2014) and Kriegel et al. (2017) at Bahir Dar, we present a new method of describing plasma depletions, which are often associated with GNSS signal scintillations. Here, we focus on describing a new plasma bubble detection algorithm and demonstrating its practical value for automatized bubble detection (and the corresponding statistical analysis). In future work, this method will be used for the automatic detection of plasma depletions and overall statistical analysis of related scintillation occurrences at GNSS measurements in the equatorial ionosphere. These stations are designed to receive and process multiple GNSS signals such as GPS (L1, L2, and L5), GLONASS (L1 and L2), Galileo (E1 and E5a), and BeiDou (B1, B2, and B3) at high data rates (e.g., Bahir Dar: 50 Hz and Tenerife: 20 Hz). The DLR GNSS receivers are connected to choke ring antennas to suppress multipath effects and allow tracking at low elevation scenarios. After removing cycle slips, key ionospheric parameters like amplitude scintillation (S 4 ), phase scintillation ( ), and slant total electron content (STEC) are calculated instantaneously and streamed together with the GNSS raw data to the central processing and controlling facility of the Experimentation and Verification network (Noack et al., 2005), operated at DLR, Neustrelitz, Germany. Sample skyplots in Figure 1 (bottom panel) show the calculated S 4 index (color coded by magnitude) for all satellites in view as a function of azimuth and elevation. As expected, enhancements of scintillation activities over the Tenerife Island GNSS station occur mostly to the south of Tenerife Island due to its position relative to the northern equatorial crest. Accordingly, the skyplots of Bahir Dar typically show enhanced S 4 values both in northward and southward direction since Bahir Dar is located in the middle of the maxima of the equatorial anomaly region (at approximately 15 • north and south of the geomagnetic equator). The occurrences of enhanced radio scintillations at both crests of the equatorial anomaly are rather typical (e.g., Hanson & Bamgboye, 1984). Nevertheless, it requires more studies to better understand the underlining physics.

Scintillations and Depletions
Strong electron density irregularities characterized by horizontal scales in the order of the first Fresnel zone may cause severe amplitude and phase fluctuations in GNSS signals. The top subpanels of Figure 2 show the observed scintillation indices (S 4 , ) at measurement station of Bahir Dar 02 (msbd02), on 16 March 2015. These observations were derived from GPS G21 and GLONASS R20 satellite signals. Near sunset, the strength of the prereversal enhancement causes equatorial spread F and plasma bubbles, in turn, are playing a dominate role in the occurrence of a strong amplitude and phase scintillations (e.g., Fejer et al., 1999;Woodman, 1970). In particular, the top and bottom subpanels of Figure 2 show a typical correlation of enhanced scintillation activities and STEC depletions observed by GPS G21 and GLONASS R20 satellite at Bahir Dar, on 16 March 2015. The depletions are believed to be caused by the Rayleigh-Taylor instability (e.g., Farley et al., 1970). It has to be noted that the used STEC measurements are not calibrated; that is, they are biased by a constant value. Thus, the analysis is not impacted by a simplifying of mapping function. Furthermore, link-related features like elevation and azimuth of considered signal paths can be used in subsequent studies to discuss the role of ray path geometry related to plasma depletion structures and associated scintillation activity. The slight time lag between enhanced scintillation and STEC depletion observed by the GLONASS R20 satellite around 21:52-22:52 LT may be derived from the interaction of GNSS link with limbs of plasma depletion, that is, attributed to the relative location and geometry of the receiver-satellite link with respect to plasma depletions in the crest of the equatorial anomaly. Generally speaking, a high correlation between the occurrence of equatorial plasma depletions and radio scintillations seems to be typical at low latitudes. Thus, more systematic and detailed studies on the relationship between equatorial plasma depletions and radio scintillations are required to interpret the occurrences and characteristics of equatorial plasma bubbles in a scientific way. The combined observations of Figures 1 and 2 indicate the concurrence characteristic of scintillations and plasma bubbles. Equatorial plasma bubbles often occur in the crest of the anomaly regions, and their monitoring can help us understand the generation and propagation of ionospheric irregularities in the low-latitude ionosphere (e.g., Valladares et al., 2001). To enable such studies by a broad database, a STEC depletion detection algorithm has been developed. This tool can be used to comprehensively demonstrate the relationship between plasma bubbles and scintillations by automatically analyzing a huge database of high-rate GNSS measurements.

Depletion Detection Method
As mentioned before, we focus on the detection of wedge-shaped structure of plasma depletions. The main descriptive parameters of such a structure are the depth and the width of the plasma depletion. The meaning of both parameters is shown in Figure 3. Hereafter, the width of the observed STEC depletion is defined as the time between the peak slopes of the wedge-shaped depletion, which means the time between the reduction phase (D) and the subsequent recovery phase (F). The STEC values marked at the well-defined position in (A), (B), and (C) are then used to define the depth of the STEC depletion, that is, the averaged STEC marked in (A) and (C) down to the minimum value marked in (B) (cf. equation (3)). Since using the observed data from a single GNSS station, a measured width of plasma bubble is characterized as a pseudowidth, which may move at an uncertain drift velocity and unknown direction. Moreover, there may be an unknown geometrical relationship between the GNSS link and the plasma bubble structure. These facts have to be considered when interpreting the results obtained by the detection algorithm proposed in this paper. To guarantee a representatively measured STEC, the observations are fitted by a polynomial function of fourth order. Its derivative can be represented by a polynomial function of third order, which will be useful to identify the three major phases (on, center, and off phases) of the plasma depletion structure. They can be expressed successively as follows: where P(t) stands for the wedged-shaped structure of plasma depletion and a 1 , a 2 , a 3 , a 4 , and a 5 , are the fitted coefficients. When computing the time derivative of the polynomial function described by equation (1), we get where P ′ (t) stands for the slope of the plasma depletion and b 1 , b 2 , b 3 , and b 4 are the coefficients of the smooth slope function. After fitting the polynomial function of fourth order to the measured STEC data, we compute the time derivative and can then define characteristic shape parameters of the plasma depletion, that is, the depth and width of the plasma depletion. These parameters can be easily deduced from STEC and time values at Points A, B, and C. In the right bottom subpanel of Figure 3, the slope values at Points D and F give a good indication about the steepness of the depletion walls. The deduced values P(t) and P ′ (t) are the key quantities for describing and characterizing the structure of equatorial plasma depletion and subsequently assist in the calculation of the depth and width of the STEC depletion. All specific positions shown in the right subpanels of Figure 3 (A, B, and C) and (D, E, and F) indicate the on, center, and off phases in the structure of plasma depletion as well as the minimum slope in the on phase, the inflection point in the center phase, and the maximum slope in the off phase in the STEC depletion structure. Referring to left subpanels of Figure 3, the peak scintillation is somewhat pronounced within edges of the STEC depletion, that is, using two points that would be located well inside the plasma depletion, rather than at "the outer lips" of the plasma depletion. This observation prompted us to develop our method for describing plasma depletions with only a few parameters that can easily be deduced from measurements obtained from an automatized data screening process within a maximum time frame of 1 hr. The amount of STEC in the on and off phases represents edges values in the STEC depletion structure. Hence, edges values can also be considered as a reference boundary between the plasma depletion and the background ionosphere. Considering the observed shape of equatorial plasma depletion, it is possible to define depth and pseudowidth of the equatorial plasma depletion accordingly. The depth (d) of the equatorial plasma depletion is the STEC difference between the edge phase and the center phase in the structure of the observed STEC depletion. It is written as follows: where STEC A , STEC B , and STEC C are the amount of STEC in the on, center, and off phases of the equatorial plasma depletion, respectively. The pseudowidth (pw) of the examined STEC depletion is defined by the time interval between Points A and C of the GNSS measurement: Taking into account the polynomial fitting procedure, it can be concretized as follows: Here, pw refers to the pseudowidth, commonly called the width of the plasma depletion, and t P ′ (t) min and t P ′ (t) max are the time of the minimum and maximum slope of the plasma depletion structure, respectively. The pseudowidth is measured in units of time indicating that it does not represent the real spatial structure of the depletion. This is due to the relative motion of the depletion and the location of the GNSS signal, which is scanning the plasma depletion depending on azimuth ( ) and elevation ( ) angle of the line of sight between the GNSS receiver and satellite. Thus, if a component of the depletion drifts vector (e.g., Kriegel et al., 2017) is in line with the velocity vector of the ionosphere piercing point, pw increases. In case, the depletion moves in the opposite direction as the piercing point moves, pw decreases. It is assumed that these motion related effects cancel out gradually by the better employment of the statistics. Analogously to scintillation and STEC depletion shown in the left subpanels of Figure 3, the top right subpanel of Figure 3 demonstrates the predominant structure of plasma depletion obtained from measurement (blue), Savitzky-Golay filter (yellow), and polynomial fitting technique (green). The bottom right subpanel of Figure 3 indicates the slopes of the STEC depletion (P ′ (t)), in particular marks the peak slopes as well as the inflection point (center) of the plasma depletion. From a technical point of view, the specific structure of plasma depletion is obtained by applying a sliding window computed every minute. The maximum length of the sliding window is around 1 hr to cover a full depletion. This method enables us to capture the required shape of STEC depletion from the overall observation (cf. right top subpanel of Figure 3). Edges and sharp slopes of the STEC depletion can be seen simultaneously at a given time axis. It is assumed that the sharp slopes of the STEC depletions can be attributed to enhanced scintillation activity in the defined time frame. In practice, the edges of the STEC depletions can be considered as a boundary between the plasma depletion and the background ionosphere. When considering the fitted polynomial function of fourth order in comparison with the observed data from the right top subpanel of Figure 3, it shows a good approximation and their root-mean-square deviation varies in the range from 1.5 to 4.4 TECU. These deviations indicate sometimes wave-like structures, which could be extracted as a by-product to characterize wave-like changes of plasma caused by other perturbation processes such as gravity waves. The simultaneous detection and analysis of traveling ionospheric disturbances (TIDs) may even help to study the relationship between gravity waves and the occurrence of plasma bubbles.

Application of the Method
The developed method has been applied to GNSS measurement data and has proven to be able to derive the structure parameters of plasma depletions. Basic requirements such as S 4 > 0.2, depth d ≥ 10 TECU, pseudowidth pw ≥ 15 min and slope P ′ (t) ≥ 10 mTECUs −1 are set up and applied on measurement data in the first study. So the influence of TIDs characterized by lower TEC amplitude can be avoided in the observation. The preliminary setting of these parameters will be proven in future studies and modified if required. The peak slopes of the STEC depletion are measured using a metric system [mTECUs −1 ], where 1 mTECU = 10 −3 TECU. To demonstrate the applicability of the developed tool for the detection of plasma depletion, we routinely analyzed additional data sets obtained from Lomé, Dakar, and Lima (see Figure 1). Based on the data availability, the percentage occurrence of plasma depletions with a depth level equal to or greater than 10 TECU amounted to 17% of days in 2015 at Bahir Dar, 7% of days in 2015 at Tenerife, 21.6% of days   Figure 4 illustrate that the maximum occurrence of plasma depletions was observed at the equatorial ionosphere during equinox months. Unfortunately, frequent power outages caused data gaps at the Bahir Dar station, mostly happened between September and December in 2015. These outages might be responsible for an incomplete observation of peak plasma depletions between September and December in 2015. Taking the reduced data availability into account, the results are in line with other observations showing maxima of plasma bubble occurrence rate during equinox months (e.g., Blanch et al., 2018;Magdaleno et al., 2017). These plasma depletions, in turn, are associated with the maximum occurrence of scintillations in the African low-latitude ionosphere, as also indicated in the studies by Hlubek et al. (2014) and Akala et al. (2015). Our method for the automatic detection of plasma bubbles is comparable to a great extent with the detection techniques developed by Nishioka et al. (2008) and Blanch et al. (2018). As described by Nishioka et al. (2008), the differential value, |R ev − R da |, was used as index of plasma bubbles activity for one day at one station. Here, R da stands for the average value of ROTI during 3 hr from the noon, and R ev is the evening ROTI from sunset to midnight. When detecting plasma bubbles by the method of Nishioka et al. (2008), larger values of |R ev − R da | ≥ 0.1 have been taken to reduce the noise level of ROTI for each GPS receiver. According to Blanch et al. (2018), SIGMA (2DTEC) ≥ 0.714 TECU, disturbance duration ≥ 10 min and depth ≥5 TECUs, was considered as a threshold or requirement for the detection of plasma bubbles. Figure 4 indicates that the results of plasma bubble detection techniques convincingly show the seasonal occurrence of plasma bubbles by using our method (left subpanels), Blanch et al. (2018) method (middle subpanels), and Nishioka et al. (2008) method (right subpanels), respectively. As can be seen in Figure 4, the number of plasma bubbles detected by Blanch et al. (2018) is greater than the number detected by our method probably due to its lower thresholds and also due to TIDs, which commonly interfere with bubble signatures. The right subpanels of Figure 4 show the monthly occurrence rate of plasma bubbles and sunspot number. Generally speaking, the detection results observed from these three methods confirm that the occurrence of plasma bubbles peak near to equinox months in 2014 and 2015 (cf. Figure 4). The histogram of TEC depletion depths shown in Figure 5 agrees quite well with findings by Rama Rao et al. (2006); their study also indicated that mostly TEC depletions with durations ranging from 5 to 25 min and magnitudes from 5 to 15 TECUs are associated with L-band scintillations. It should be mentioned that in addition to the increased scintillation activity also the positioning accuracy of GNSS is affected by several meters. Compared with other methods referenced here, our detection method that describes depth and width of plasma depletions and has the advantage that it relies on smoothed STEC data, which are not impacted by mapping function errors. Our method can give the benefits for ionospheric scintillation studies when the derived TEC depletion parameters are available (see example in Figure 5). Furthermore, our method clearly indicates when the satellite-receiver links enter and leave a depletion zone by the derived slopes (Points A, C and D, F, respectively). This enables us to perform further studies on the relationship between ionospheric plasma bubbles and radio scintillations and their impacts on GNSS signals.

Conclusions
Simultaneous observations of enhanced scintillation activities and STEC depletions have shown their correlation across the GNSS links. Based on these preliminary observations, a new method has been developed that is able to describe the main features of the structure of equatorial plasma depletions, including a characterization of the depth and the width of the depletion structure. We found that depth, width, and peak slopes of the STEC depletions may range from 10 to 40 TECU, from 15 to 68 min, and from −10 to 40 mTECUs −1 , respectively. The presented results have been validated and compared with other GNSS plasma depletion detection techniques indicating a qualitative agreement. The deduced structure parameters such as depth, width, and slope of STEC depletions can help us to study the 3-D shape of plasma bubbles and the relationship between plasma depletions and scintillation activities. Having a better understanding of this relationship will be helpful in nowcasting and forecasting ionospheric scintillation activities of GNSS signals. The proposed method will be applied in more comprehensive studies for automatic detection of plasma depletions and related radio scintillation occurrences in the equatorial ionosphere.