Automated Detection and Tracking of Equatorial Plasma Bubbles Utilizing Global‐Scale Observations of the Limb and Disk (GOLD) 135.6 nm Data

Equatorial plasma bubbles (EPBs) are regions of depleted plasma within the ionosphere that form during post‐sunset hours near the magnetic equator. EPBs tend to align with local geomagnetic field lines, extend hundreds of kilometers along geomagnetic longitudes, and thousands of kilometers along geomagnetic latitudes. EPBs can attenuate lower frequency radio waves, such as HF radar. Smaller features, usually 10s of meters, along the EPB's wall can interfere with centimeter scale wavelengths, such as GPS, via Bragg scattering. Statistical analysis of EPBs can further understanding of their occurrence and behavior. The current study utilizes Global‐Scale Observations of the Limb and Disk (GOLD) 135.6 nm nightglow data from 5 October 2018 to 30 September 2022. GOLD has a unique perspective from geostationary orbit, allowing consistent viewing of nightglow and structures over the Americas and Atlantic. An EPB detection method is developed and used to generate a database of occurrences. Occurrences are used to calculate EPB drift speeds and separations. All data is looked at as a whole, including any possible solar or geomagnetic disturbance periods. Seasonality in occurrence rate is evident. Occurrences peak during December solstice and minimize during June solstice for longitudes seen by GOLD. Within GOLD's view, higher occurrences are seen to the west during December solstice and to the east during June solstice. EPB drift speeds and separations show consistent distributions regardless of magnetic latitude region, geographic region, or season. This suggests EPBs behave consistently and regularly once formed, at least on spatial scales observed by GOLD.

EPBs tend to align with local geomagnetic field lines.From a geomagnetic reference frame, EPBs will range from several to hundreds of kilometers along magnetic degrees longitude and thousands of kilometers along magnetic degrees latitude (Immel et al., 2003).These structures also rise and stretch from the lower F-region of the ionosphere to well over 1,000 km (Comberiate & Paxton, 2010).EPB drifts are generally west to east with the background zonal plasma drift (England & Immel, 2012;Fejer et al., 1991;G. Liu et al., 2013;Makela, 2006).These regions of depleted plasma tend to have a decreasing velocity as we trace them upward in altitude which causes a perceived westward tilt (Aa et al., 2020).Variations in drift speeds along magnetic latitudes can lead to EPBs taking a C or reversed C-shape (Kelley et al., 2003;Kil et al., 2009;G. Liu et al., 2013;Mendillo & Baumgardner, 1982;Mendillo & Tyler, 1983;Woodman & La Hoz, 1976); however, this has also been attributed to neutral zonal winds at F-region altitudes (Anderson & Mendillo, 1983;Martinis et al., 2003).The observations required to develop an understanding of EPBs come from remote sensing and in situ devices.Ground-based airglow and radar provide bottom side information, in situ measurements provide data both within the EPBs and of the surrounding plasma, and global GPS TEC data allow continuous viewing of EPBs.
Large scale statistical analyses over many years consisting of thousands of EPBs have been done previously and provide an understanding of the seasonality, frequency, and general behavior of EPBs during various periods of the solar cycle.One of the largest of these studies by Sobral et al. (2002) showed seasonal variations and increases in EPB occurrence rates with solar activity from 1977 to 1998 over Cachoeira Paulista (45° west, 33° south, and 28° southern dip).
A study from Gentile et al. (2006) used plasma density data from DMSP to detect EPBs from the surrounding plasma and define EPB occurrence rates with respect to total orbits from 1989 to 2004.Huang and Roddy (2016) used C/NOFS ion density data from 2008 to 2014 to detect EPBs.The lower inclination orbit of C/NOFS ensured minimal variations in latitude from orbit to orbit which allowed for the tracking of EPBs and the ability to calculate EPB drift speeds and to see how they were affected by solar and geomagnetic activity.Finally, Wan et al. (2018) used Swarm electron density to detect and analyze the generation of EPBs.The near-polar orbit of Swarm allowed for EPB characterization and analysis of formation along the latitudinal direction.
The current study further develops this understanding using Global-Scale Observations of the Limb and Disk (GOLD) 135.6 nm airglow nightside data by taking advantage of GOLD's location in geostationary orbit (GEO) to detect and track EPBs consistently, after their formation, throughout the night.A systematic and automated method is developed to process GOLD data and detect EPBs therein.The automated process lends itself to the processing of many years worth of data with no external input while the systematic nature of the method ensure consistent results within a given set of data.The method is also developed and explained in a way that allows it to be altered and changed to fit the data and requirements of anyone who may choose to implement it.A detrending scheme developed by Pradipta et al. (2015) is used to pick out and generate a database of EPB occurrences.This detrending scheme minimizes the stretching and skewing of data that can occur from other detrending or averaging techniques.From there, a process is defined and used to identify EPBs from one GOLD image scan to the next for drift speeds to be calculated.Unique EPBs are defined as EPBs that appear in multiple scans and are deemed the same EPB by way of drift speed calculations and are counted.This allows for the detection and tracking of multiple EPBs per scan across dozens of scans per night over the same geographic region over years.These EPBs are also used to calculate interbubble separations.EPB occurrences are analyzed by season and geographic regions.Analysis of EPB drift speeds and separations by geographic region, magnetic region, and season provides a view of EPB behavior within GOLD's view.

GOLD Mission and Data
The GOLD mission has a unique perspective on Earth's TI from its vantage point in GEO above −47.5°longitude.This provides GOLD with consistent data on the spatial and temporal evolution of the dayside and nightside TI over the American, Atlantic, and west African regions simultaneously via routine observations of the Earth limb and disk (Eastes et al., 2017(Eastes et al., , 2019)).The GOLD imager consists of two identical and independent channels that can operate individually or simultaneously to collect 132-162 nm emissions from the TI by scanning across the disk and/or limb at regular intervals as little as 15-30 min, dependent on operations.Information on the temperature and composition of the TI system can be gleaned from these measurements and allow for the separation of temporal and spatial variations due to the consistent nature of the observations (Eastes et al., 2017).Nightside data is used to observe variations in radiance given off by radiative recombination of   + with electrons.These variations correspond to structures in the low magnetic latitudes such as the equatorial ionization anomaly (EIA) and EPBs (Eastes et al., 2019(Eastes et al., , 2020)).Nightside scans start east of the dayside terminator along the eastern edge of the Earth disk and progress westward following the terminator with alternating north and south scans with a temporal cadence of ∼15 min.Scans are taken simultaneously with both channels later in the evening.Each pixel within a scan has a spatial resolution of ∼95 km at nadir, ∼112 km nearer the limb, and are ill-defined at the limb (Eastes et al., 2017(Eastes et al., , 2020;;McClintock, Eastes, Beland, et al., 2020;McClintock, Eastes, Hoskins, et al., 2020).
The automated EPB detection and tracking algorithm is implemented on data from 5 October 2018 to 30 September 2022 for statistical significance.Nightside (NI1) level 1C (L1C) version 4 data is used which contains geographic coordinates for corrected photon particle counts and calibrated radiance values.This data is collected east of the dayside solar terminator and contains local times in the range of 19:00-22:00 local time.Scan cadence has varied throughout the mission with the advent of special operations and changing mission parameters.The method has therefore been designed to operate effectively with varied frequency of data.

Method Outline
The method for identifying the locations and characteristics of EPBs in GOLD images has the following major steps: 1. Data are interpolated onto a regular grid in quasi-dipole (QD) coordinates to identify and later track features about the magnetic equator (Sections 3.2-3.4). 2. Individual north and south scans are combined into full images and the region of interest (ROI) is isolated (Sections 3.5 and 3.6).3. EPBs are isolated from the background ionosphere and their properties including location, width, and intensity are recorded (Sections 3.7-3.11).

Cubic Spline Interpolation
The airglow data from GOLD is in a coarse irregular spatial grid in geographic coordinates, but the method requires data to be in a regular grid in QD coordinates.This necessitates a two dimensional interpolation in array space for the Nearest neighbor interpolation (NNI) that follows.A simple two dimensional clamped cubic spline is implemented to accomplish this.The top row of Figure 1 shows the location of the original data points and the new interpolants in red and black, respectively, along with a zoomed in view in the middle row of Figure 1.
The number of interpolants between data points in array space used is defined by a scale value.Increasing the scale value can increase compute time substantially.A scale value of 10 in used in this study as higher values lead to substantially higher compute times with little gain in accuracy or precision.This potentially adds some artificial variation from data point to data point, but when compared to the scale of the EPBs of interest, this was deemed a sufficient method to fill in gaps in the NNI that follows.

Nearest Neighbor Interpolation
NNI is implemented for its simplicity and robust nature.This allows for simple transformation to a regular grid for analysis that comes later.NNI is done in half degree boxes with pixel centers at every half degree in longitude and latitude.The boundaries of NNI are shown in the bottom row of Figure 1.Note the regions that would otherwise be void of data without the interpolants from the cubic spline.This NNI algorithm is implemented on each GOLD scan individually.Results of NNI are shown in Figure 2a.

Quasi-Dipole Coordinate Conversion
EPBs tend to follow the EIA and magnetic field lines.The EIA is a region of increased plasma density that forms about the magnetic equator.This phenomenon can be seen in Figures 2a and 2b in geographic coordinates and 10.1029/2023EA002935 4 of 27 magnetic coordinates, respectively.The EIA appears as two semi-parallel bright regions with positive slope from west to east in Figure 2a and with minimal slope in Figure 2b.As such, conversion to magnetic coordinates will significantly simplify the detection algorithm and analysis on the back-end.The magnetic coordinates used in this analysis are quasi-dipole (QD) coordinates, similar to that used by Karan et al. (2020), and are implemented with half QD degree resolution north and south.QD coordinates with a reference altitude of 300 km are selected because 135.6 nm airglow peaks near this region (Laundal & Richmond, 2017).Figure 2b shows the results of QD conversion.

Image Pairing
The GOLD imager consists of two identical and independent channels: A and B (Eastes et al., 2017).Either one or both channels are used for nightside Earth disk observations.When one channel is used, it will alternate scans of the northern and southern hemispheres with each integration period (Eastes et al., 2017).Conversely, both north and south scans can be done simultaneously when both channels are available.This requires forming images from pairs of scans to generate a continuous time series of the nightside ionosphere.
An example of the NNI, QD conversion, and data pairing for 28 November 2018 at 23:10 UTC are shown in Figures 2a and 2b.Note how the EIA follows the magnetic equator northward and eastward in the region of the South Atlantic Anomaly (SAA) leading to a positive slope (latitude/longitude) from west to east and how this is translated to a slope of approximately 0 (QDlatitude/QDlongitude) in QD coordinates.GOLD image scans are sorted by UTC within a given night of data.This ensures images are formed from subsequent scan times.The first scan of the night is referred to as t 0 , the second referred to as t 1 , up to t i−1 for the last scan of the night where i is the number of scans for a given night.The first scan pair will then have time steps t 0 and t 1 , the second uses t 1 and t 2 , up to the final pair of t i−2 and t i−1 .This will generate up to i − 1 scan pairs as subsequent north-north or south-south pairings are ignored.When looking at scan pairs, the earlier of the time steps will be referred to as t 0 and the later as t 1 .This is for convenience and to acknowledge both scans are at different times and scan regions.
Checks are implemented during this process to ensure the images formed are proper North/South pairs with no significant gaps in time or data.An image is skipped if two scans have the bulk of their data along similar geographic latitudes.This prevents the pairing of North/North images and South/South images.A similar check implemented for geographic longitudes prevents scans that are subsequent in time but with the bulk of their data significantly distant from each other in longitude from being paired.Aside from a lack of southern scan data during early 2019 due to an issue with channel B, these issues are infrequent in the bulk of data.

Isolating the Region of Interest
A minimum solar zenith angle (SZA) is implemented to prevent any dayglow and early evening data effected by the solar terminator from creeping into the data analyzed.This ensures only steady state nightside data is utilized.The minimum SZA cutoff value is calculated using the following equation, where  max is the maximum SZA in degrees, A is altitude, and    is Earth's radius with the later two in kilometers.It can be shown through spherical trigonometry that  max = 90 • on the surface of the Earth (A = 0) and that an increase in altitude increases  max .Note that 1 assumes ideal parallel light rays through all possible layers of the atmosphere and ionosphere and ignores refraction or scattering.Dayglow from oxygen (135.6 nm) is brightest below 200 km (Meier, 1991) 2c contains a bright spot in the southwest region due to the solar terminator.This bright spot is removed with the implementation of the SZA cutoff shown in red.
Occasional polar effects, such as aurora, can be viewed through nightside disk scans.QD cutoffs of ±22 QDdeg latitude are implemented to isolate the low latitude ROI for EPB detection.The SZA and QD cutoffs are shown in Figures 2c and 2d.Note the two bright regions in the north-west and south-west extents of the geographic image of Figure 2c that are removed in the QD conversion.These bright regions could potentially lead to false EIA tracking and false EPB detections.
EIA crests that go beyond the ±22 QDdeg latitude bounds will lead to EPBs that are missed while EPBs that "slice" through the EIA will leave a larger impact to difference data calculated later for a more confident detection.The rise in solar activity in 2022 and beyond could also lead to some missed detections with shifting EIA crests.This was deemed sufficient given the time period of data being analyzed and that the EIA crests lie within these bounds for the vast majority of data looked at in this study.

Rolling Barrel Average for Detrending
GOLD sees EPBs as dips in the brightness of 135.6 nm airglow which can be seen in Figures 2a-2d (Eastes et al., 2017).These dips in brightness are most clearly seen within the EIA (Karan et al., 2020).A detrending process is used to highlight dips in brightness by way of a difference image which will correlate with the EPB locations.The averaging method used for detrending is referred to as a rolling barrel average (RBAVG) because it simulates a barrel of some radius rolling over data like a surface.This technique was developed for the detection of ionospheric anomalies by Pradipta et al. (2015) and implemented on TEC data.This study applies it to non-TEC data for the first time.This detrending method was chosen to maintain the integrity of the data being detrended without the need for further corrections: local or otherwise.Common detrending techniques, such as simple smoothing, stretch or skew structures of interest and lead to potential false negatives.The premise of the RBAVG is a barrel with sufficiently large radius will necessarily "roll over" dips in data, which can correlate to structures of interest.Specifics of the geometry are shown in Figure 3.
Following Pradipta et al. (2015), at any given data point, RBAVG looks ahead to all data points along the ROI and calculates the δ for every point therein.This δ is the angular distance between the points on the barrel that would roll over a given data point within the ROI (see Figure 3).The data point with the minimum δ is taken as the next pivot point (PP).This process assumes an ordered, not necessarily regular, domain.The data point ahead of the current PP is taken as the next PP if no data points are within the ROI which contrasts from Pradipta et al. (2015) taking the closest point outside of the ROI as the next PP.This is done to prevent skipping of data within an ROI due to noise or coarseness of input data.Once all PPs have been defined, a simple linear interpolation is implemented from one PP to the next.The curve generated is the detrended data.
Figures 4a and 4b show a smoothing average and RBAVG, respectively, implemented on a simple sine function.The smoothing average, window of 1.5π, requires a shift by a constant to match the data.More complex data with varying peaks and valleys would necessitate more localized corrections.The RBAVG, with R 0 = 1.5π as defined in Figure 3, follows the sine curve perfectly while skipping over gaps.
The detrended data are differenced with the input data to generate the red curves in Figure 4.Note how the peaks and gradients of the difference data in regions of depletions are more pronounced with the RBAVG example compared to the smoothing example.Figure 4c shows an example of RBAVG implemented on QD data from Figure 2d along constant QDlatitude of 10 QDdeg with R 0 = 4 QDdeg.This value for R 0 is used to ensure it is larger than most widths of EPBs expected to be detected.Note the three largest peaks of the difference image at 27, 32, and 44.5 QDdeg longitude and the

Equatorial Ionization Anomaly Crest Detection
EPBs are more readily seen and detected along the EIA due to the contrasting brightness of EPBs' dim depleted plasma and the brighter surrounding plasma.These sharp contrasts provide a means to detect the EPBs through the optical data.To take advantage of this, EIA tracking is implemented for the detection of EPBs to the north and south of the magnetic equator.A combination of a running average and Gaussian fit is implemented along the QD northern region and QD southern region independently to accomplish this.Data is split from the magnetic equator to the northern and southern QDlatitude cutoffs defined previously and looked at individually.EIA  tracking takes advantage of the fact that a meridional cross-section of the EIA generally appears Gaussian-like in the magnetic north and south.A 5.5 QDdeg wide window from QD west to QD east is centered about each QDlongitude.An average is then taken along QD latitudes within this window to smooth the data along QD latitudes.An example of these windows in the QD northern and southern regions centered about a QDlongitude of 30 QDdeg and the subsequent average can be seen overlaid in orange in Figure 5c.This width in QDlongitude is comparable to the widths of EPBs and would ensure the EIA tracking is not heavily skewed by their presence.Next, a Gaussian fit is implemented along the QDlongitude the average was taken.This generates a simple curve from which a max value can be taken.An example of this Gaussian fit taken about a QDlongitude of 30 QDdeg in both the QD north and QD south can be seen in Figure 5d.The peak along the northern Gaussian is taken as the northern EIA crest peak and similarly for the southern Gaussian.Implementing this across all QDlongitudes results in EIA crest peak locations that can be used as a starting point for EPB detection.Figures 5e and 5f show the EIA tracked in this manner in orange over QD data and difference data, respectively.EIA detection and difference data are the basis for EPB detection.

Plasma Bubble Detection
A single row of data is not sufficient to accurately determine what is a EPB in a given region due to the signal to noise level (see Figure 4c).Implementing summations along QDlatitudes can highlight bright spots within difference data, which correspond to features aligned with the geomagnetic field.Difference image sums are taken over the EIA peak and 4.5 QDdeg magnetically equatorward.This is done to more accurately pick out real detections, but also to ensure as many EPBs that have not quite reached the EIA crest are accounted for; however, this does not capture all such cases.Additionally, summation along the magnetic equator are taken from −4 to +4 QDdeg latitude.This yields three regions of EPB locations: northern EIA, southern EIA, and magnetic equator.The increased range for equatorial summation is to compensate for generally less intense nightglow along the magnetic equator.The QD northern, equatorial, and southern summation regions and the resulting southern sums are shown on Figures 6a and 6b.
The sums are run through a 2.5 QDdeg wide triangular convolution filter.Filtering in this manner smooths regions around shallow peaks and prevents double counting EPBs.A simple peak detection algorithm then picks out potential EPB locations.Each point is compared to all points 1.5 QDdeg west and east of itself.If it is larger than all other values, it is a peak, else it is not.Similarly, all minimum values between peaks are picked out as valleys.The results of the peak detection on both the raw sums and filtered sums for the southern region are shown on Figure 6b.
Uncertainty analysis provides the final protection against false positive detections.Both systematic and random uncertainties from the L1C data files have been propagated throughout the process.The summation in quadrature of these uncertainties is used as the measure of uncertainty: σ.Peak uncertainties are compared to the closest valley uncertainties to the east and west.If a peak minus its 1.5σ is less than both valleys plus their respective 1.5σ values, then the peak is considered too shallow and not a EPB detection.If a peak lies on the edge of the data, it is compared to the closest valley.1.5σ is used as this is effectively a 3σ uncertainty filter.Note the EPB detections and rejections on Figure 6c.The western-most and two detections between 35 and 40 QDdeg longitude have been deemed too uncertain to be EPBs when compared to the relevant valleys.Two detections at 26.5 and 44.5 QDdeg longitude are well within uncertainty bounds on either side.A detection at 31 QDdeg longitude fails the uncertainty check to the west, but passes to the east and is kept.The eastern-most detection only has a valley to the west of itself and passes its single uncertainty check.
The NNI data lead to EPBs that are detected on 0.5 QDlongitude grid.To identify the location of the center of each EPB, Gaussian fits are applied to each detection using QDlongitudes 1.5° to the east and west.The peak location of the resulting Gaussian is taken as the final QDlongitude for a given detection (Figure 7a).Note how each Gaussian lies within the domain of its respective detected peak.The final results of the detection algorithm in geographic coordinates and QD coordinates are shown in Figures 7b and 7c, respectively.EIA crest locations are used as the northern and southern EPB QDlatitudes and all EPBs detected along the magnetic equator are assumed to be at 0 QDdeg latitude for simplicity.

Determination of Plasma Bubble Width
The last thing determined before generating the final database is a measure of EPB widths.The first location to the magnetic east and west of each EPB lesser or equal 75% of the detected peak filtered sum are defined as the full width at 75% max.These are generated such that the eastern QDlongitude minus the western QDlongitude provides a measure of EPB width that can be quantified and analyzed.The value of 75% is chosen as anything less than 75% leads to EPB width extents that overlap with neighboring EPBs whereas greater than 75% leads to ill-defined and false reported widths.
Note that these measures of width are most ideal when EPBs follow lines of constant magnetic longitudes due to the summation along QDlatitudes (Section 3.9).The widths will be overestimated should the EPBs take on a C-shape or reversed C-shape as shown in Mendillo and Tyler (1983) and G. Liu et al. (2013).

Database Generation
The following values are saved for a given image pair for each EPB once all bubbles have been classified for each magnetic region; Gaussian coefficients for EPB center location, QDlongitude and latitude, year, day of year (DoY), UTC for both scans that comprise an image, intensity sum and filter values, random and systematic uncertainties pre and post convolution filter, western and eastern 75% max measures of width, and each individual image's geographic scan region: north or south, SZA, and the intensity of crest locations.The full 4 years is available, see open research statement.

Selection, Drifts, Counting, and Separation
Two things prevent accurate drift calculations of EPBs.Firstly, GOLD images can share a north or south scan from one image to the next leading to a single double count at best: multiple double counts at worst.Second, there is no definitive way to show a given EPB is the same one image to the next.The first is simple to check for.Various approaches have been used for the second and are dependent on the data set.S. H. Park et al. (2007) took advantage of IMAGE-FUV's low temporal cadence compared to the motion of EPBs (approximately 2 min) to detect EPB locations and drift velocities by looking at keograms to develop their automated EPB detection algorithm and calculate EPB drift speeds with low spatial resolution: ∼100 km at nadir at apoapsis (Mende et al., 2000).Huang and Roddy (2016) used 6 years of high resolution in-situ data from C/ NOFS with a temporal cadence dependent on its orbital period of ∼100 min to identify individual substructures within each EPB and calculate drifts.GOLD's nightside data lies somewhere in the middle.From Section 2, the nightside scans have a temporal cadence of ∼20 min with more frequent scans in recent data products and a spatial resolution of ∼95 km at nadir (McClintock, Eastes, Hoskins, et al., 2020).This temporal cadence allows for consistent and more frequent tracking of bubbles from one scan to the next with less uncertainty within the same night and does so from a fixed vantage point.Karan et al. (2020) showed how GOLD's nightside data could be used to detect, track, and calculate drifts of EPBs by manually selecting and analyzing dips in radiance which correlate with regions of depleted plasma.S. H. Park et al. (2007) and Huang and Roddy (2016) showed the power in implementing an automated method for detection, analysis, and statistics regardless of data.
The automated process implemented for this research and the relevant sections are as follows: 1. EPBs are selected from distinct images (Section 4.1).2. Drifts are calculated from one image to the next (Section 4.2). 3. EPBs are counted (Section 4.3).4. Double counts removed (Section 4.4). 5. EPB separations are calculated (Section 4.5).6. Uncertainty calculations used (Section 4.6).

Plasma Bubble Selection
Comparing one image's EPBs with the next image's EPBs must be done to calculate EPB drifts and counts.Care must be taken not to compare the same scan with itself from one image to the next as can be seen in Figures 8a-8d and must be done in such a way that accounts for changing scan cadence of both channels and special operations over the course of the GOLD mission.All southern region EPBs would be compared with themselves if Figure 8a were compared with Figure 8b as both have the same southern scan.Comparing Figure 8a with Figure 8c would ensure unique scans in both images.The same logic can be applied if the current image were instead Figure 8b.This is accomplished by comparing the UTC and geographic region of both scans that comprise an image in a given night.If there exists a common UTC and north/south scan, then the next image is checked until a valid image comparison can be done.For example, if Figure 8a is the current image, then it will be first compared to Figure 8b.Both Figures 8a and 8b have a common UTC south scan, and so Figure 8a is then compared to Figure 8c.Both have unique scans from each other and so EPBs can be used for later drift and tracking calculations.The same will then be done for Figure 8c and subsequent images.The same can be done if Figure 8b happened to be the starting image.Figures 8b and 8c have a UTC north scan in common, and so Figure 8b would be compared with Figure 8d, and so on.This process is accomplished more rigorously as follows.

Drift Calculation
EPBs are tracked one at a time by comparing a single EPB in one image with all EPBs in the next image that lie slightly west and fully east of its last known position.This is done for each magnetic region individually.This method of tracking was chosen over a full image correlation analysis method for several reasons.EPBs tend to drift magnetically east while GOLD scans move progressively westward with UTC with an integration period of ∼15 min.This leads to image comparison time differences of ∼15 min when both channels are actively scanning, such as those shown in Figures 8a-8d.When only a single channel is actively scanning, then image comparison time differences are typically ∼30 min.The cadence of scans that form the images also changes over the course of the GOLD mission for changing mission operations and special operations which can lead to more inconsistency with time differences from image to image can lead to image comparison time differences of up to ∼60 min.
Features of interest are able to propagate, change shape, and disappear from view within this time.The large and varying time steps along with non-ideal overlap in longitude from image to image therefore does not lend itself to a full image correlation analysis which assumes drifts are largely consistent and features persist in a regular manner.
All EPBs are initialized with a drift of zero.Drifts are calculated by comparing EPBs from   and  +1 in a QDwest to QDeast manner by magnetic region.This is done to take advantage of the generally eastward drift of nightside EPBs (England & Immel, 2012;Huang & Roddy, 2016) The final drift magnitude is simply, where v j,n+1 is positive eastward (negative westward) as determined by 3.
The highest previously noted drift speed is ∼190 ms −1 by Mendillo and Baumgardner (1982) and Sinha and Raizada (2000).When converted to degrees per hour at higher latitudes, this can range from 7 to 7.An example of these calculations is provided in Figure 8e  It is important to note these drift calculations are taken from the EIA peaks at   's to  +1 's location as defined in Section 3.9.Therefore, this method does not provide information on gradients in EPB velocity with respect to QDlatitude.

Plasma Bubble Counts
This method requires EPBs to be seen in at least two images to be counted and tracked.This ensures only EPBs that persist longer than a single image integration period are counted.All EPBs are initialized as having a count number, c, of zero (uncounted) with, where "count" is an index value for the database accessing EPB count numbers and c 0 = 0 is the initialized EPB count value.EPBs are counted in order drift calculation as defined in Section 4.2.
If 7 is false, then a check on  (count) and  +1(count) is done and c values are updated as follows, where  max is the current max c value in the database thus far.

Double Count and Untracked EPB Removal
EPBs may be detected several times across several images and these individual detections are heavily biased by GOLD's scan cadence: more/fewer scans lead to more/fewer EPB detections.This does not necessarily imply more/fewer unique EPBs.The removal of double counts and untracked EPBs ensure only unique EPBs viewed in at least two subsequent images are kept.Section 4.3 ensures the removal of all double counts and untracked EPB occurrences is now trivial and occurs following drift calculations and count assignments.With all EPBs labeled with count values greater or equal one, simply remove all EPBs in the database that are still labeled as c 0 .Implementing this removal ensures only tracked EPBs are kept for a given night.But will only be counted as a single EPB.This ensure an increase/decrease in scan cadence, and therefore an increase/decrease in total detections, does not heavily bias total EPB counts.

Plasma Bubble Separation
Similar to prior sections, all EPB separations are initialized as zero.Plasma bubble separation is calculated in QD and geographic coordinates by looking through all EPBs in a single image as selected in Section 4.1 by magnetic region.It is worth noting that some separation values may be overestimated due to the removal of initial detections (Section 3.9) and EPBs with no associated drift speed (Sections 4.3 and 4.4).All results will be interpreted with this in mind.
This can be seen in Figure 8f using EPBs from Figure 8c.The southern EPB at ∼15.5 QDdeg QDlongitude has a separation calculated with the EPB at ∼20 QDdeg QDlongitude.The separation is ∼600 km and is plotted at the western EPBs location as ∼15 QDdeg QDlongitude and ∼600 km.This is then done for the EPB at ∼20 QDdeg QDlongitude and ∼24 QDdeg QDlongitude and so on.The eastern-most EPB along the southern EIA is shown as having a EPB separation of 0 km to show it is the eastern-most EPB in this region that has been detected and tracked.

Uncertainty Calculations
QD position uncertainties are calculated during Section 3.9 and take summation uncertainties into account.QD coordinate uncertainties are added in quadrature from one EPB to the next for both drift and separation calculations as follows, for v j,n+1 uncertainty and, values are relevant coordinate uncertainties in QDdeg and e is a conversion scalar with units    .The first term of e is the maximum nearest neighbor distance per its distance in degrees within the geographic ROI.The second term of e is the minimum QD degree difference per geographic grid spacing during coordinate transformation.e is approximately  80   for the current work.

Results
The method developed in Sections 3 and 4 is implemented on all the data from Section 2 to create a single database of EPBs.Generating a database grants ease of access and filtering for deeper analysis.
The results that follow include all data with the rise in solar activity since 2018 and geomagnetically active periods.The geomagnetically active periods that are present in the current database are deemed few enough in number that an overall statistical analysis is sufficient for the current study.This is done to develop an overall statistical view on the entire data set.Future work will engage with solar and geomagnetic active periods along with coincidences with other missions.Results are consolidated and shown as follows: 4. EPB Drift Speeds (Section 5.4), 5. EPB Widths and Separations (Section 5.5).

Plasma Bubble Detections
A total of 26,714 detections along the northern EIA, 29,537 detections about the magnetic equator, and 24,913 detections along the southern EIA are deemed as valid and used for further analysis.Of these detections, there are 7,443 tracked EPBs along the northern EIA, 8,196 tracked EPBs along the magnetic equator, and 7,136 tracked EPBs along the southern EIA.The magnetic equator is generally dimmer and leads to the possibility of more false positives than the northern and southern EIA.This is likely due to the uncertainty calculations for the L1C radiance values used to detect the EPBs: details in Section 3.1 of McClintock, Eastes, Beland, et al. (2020).
One term the uncertainties are proportional to is the total detector counts accumulated over an integration period (McClintock, Eastes, Beland, et al., 2020).This leads to peak and valley comparisons in Section 3.9 that are more shallow with smaller uncertainties.It's then possible for small fluctuations in radiance along the magnetic equator to be detected and reported as a potential EPB.Keeping only detections that lead to an EPB drift calculation has ensured that not all such radiance fluctuations are reported as EPBs, but some will inevitably persist.The aforementioned drift calculations coupled with the uncertainty checks from Section 3.9 ensure there is a measure of consistency and provides some confidence that these extra detections are valid.There are also periods of time where one or both EIA crests are visible within the magnetic equator's ±4 QDdeg QDlatitude summation region.These detections will be detected along the magnetic equator leading to a higher number of detections about the magnetic equator.
Heat maps of individual EPB detections are shown in Figures 9a-9d.Total detections are highest during December solstice (November-February), lower during the September equinox (September and October) and March equinox (March and April), and lowest during June solstice (May-August).This is corroborated by the time series plots in Figures 10a-10f.Detections during the December solstice and equinoxes are largely symmetric about northern and southern EIAs, but there exists an asymmetry within June solstice, especially over west Africa.This phenomenon is also visible in prior studies, such as J. Park et al. (2022).The phenomenon has been attributed to the sun's apparent location during the day and the solar terminator being nearly parallel to the QD field lines about the SAA.
The spikes in total detections during the 2022 solstices and equinoxes can be attributed to GOLD's change in mission operations which increased scan cadence from ∼24 scans per night to ∼36 per night.This increase in scan cadence leads to an increase in useable images for EPB detections and therefore an increase in total detections.This is accounted for by counting unique EPBs.

Seasonal Variation of Plasma Bubble Counts
After accounting for unique EPBs (Section 4.4), the spike in detections during the early months of 2022 are no longer evident in Figures 10d-10f.This shows the tracking and EPB counts are not affected by GOLD's change in operation.Using unique EPBs, GOLD sees a quasiperiodic behavior over the year and a slight general increase from 2018-2022 (Figures 10d-10f).Comberiate and Paxton (2010) and Magdaleno et al. (2017) show consistent seasonal variations with GUVI data and GPS data, respectively, in similar geographic regions comparative to the time series plots of unique EPB counts shown in Figures 10d-10f.Further, the EPB detections and counts exhibit quasiperiodic behavior.Rise of EPB occurrence rate from 2019 to 2022 are most likely due to rising solar activity over these years.Comberiate and Paxton (2010) show decreases in EPB occurrences from 2002 to 2007 while Magdaleno et al. (2017) show a similar solar dependence; both corroborate these results.
It is worth noting the current research is using unique EPBs per night via tracking of EPBs to develop these statistics.In this case, the more unique EPBs detected imply a more active night and allows for a night by night time series to be generated.Later sections also look at the percentage of EPBs over an entire subset of the population split by season and region.Prior researches have used different definitions to develop their statistics.For example, Martinis et al. (2021) defines occurrence rate as number of nights at least one EPB was detected per days in the month of observation.In this case, 4 EPBs in a single night of a 30 night month would still count as 1/30 of a given month's occurrence rate.The varying definitions of occurrence rates to develop statistics can lead one to defining things like active months or periods differently.

Plasma Bubble Characteristics by Geographic Region
Geographic regions are defined to identify any variation in EPB characteristics by location as follows: South American (SA), Atlantic (AT), and west African (AF) sectors from −90° to −34°, −34° to −18°, and −18° to 15° longitude, respectively.The eastern and western limits are chosen arbitrarily outside of GOLD's view.The aforementioned regions were chosen as the geomagnetic equator has a sharp positive slope from east to west with QDlatitude over SA, shallower and decreasing slope over AT, and nearly zero slope over AF.Figures 11a and 11b show EPB detections across the whole database and by season, respectively, by geographic region.More EPB detections are visible to the west in Figure 11a, but Figure 11b shows this largely depends on season.Martinis et al. (2021) looked at EPB occurrence and formation with GOLD data during November 2018 and March 2019 in several geographic regions.They note higher occurrence rates in their West-American and Brazil-West Atlantic sectors, which coincide with the defined SA region, and lower occurrence rates in their East Atlantic-West-African sector, which coincide with the defined AF region.The December solstice heat map shows an overall agreement with Martinis et al. ( 2021)'s results over the entire database; however, the spring equinox     Huang and Roddy (2016) reported average drift speeds around ∼125 m s −1 for the equinoxes and December solstice when F10.7 was 135sfu and just under 100 m s −1 when F10.7 was 115sfu during the local times that GOLD's nightside data is collected.Huang and Roddy (2016) also showed a dependence of drift speed on solar flux, namely a linear increase in EPB drift speeds as F10.7 increases.Average F10.7 for the GOLD data being analyzed was ∼85sfu with average drift speeds ranging from 70 to 110 m s −1 .
Westward drifts mentioned in Taylor et al. (1997) are few in number, but persist for upwards of 90 min.While these are well within GOLD's temporal resolution, these occur during and after local midnight and during low and high K p .Most reported westward drifts of EPBs during the local times GOLD observes tend to appear during active periods such as those reported by Makela (2006).J. Park et al. (2021) used ICON nightside data to show how plasma within EPBs correlates with the surrounding plasma.While the EPBs tend to have a more westward drift than the surrounding plasma along the mid to low latitude regions, there does not appear to be any evidence of EPBs drifting westward.This shows that westward EPB drifts are possible and within GOLD's ability to view.However, GOLD is not able to measure the F-region winds.Therefore, changes in solar activity, geomagnetic activity, or anomalous F-region winds are all possible causes of the westward EPB drifts seen.
As a final note, drift speeds in Figures 11c-11h and Figures 12a-12d show a number of drift speeds greater than 200 m s −1 .The highest previously noted drift speeds used in Section 4.2 as a maximum cut off were ∼190 m s −1 by Mendillo and Baumgardner (1982) and Sinha and Raizada (2000), but when converted to degrees per hour at higher latitudes allow some higher drift speeds to pass at lower latitudes.Aside from June solstices when data is more sparse, these generally account for around %6-%7 of calculated drift speeds and are recognized as outside of the norm from previous research.

Plasma Bubble Widths and Separations
Widths in the range of 100-400 km are detected and displayed in Figures 13a and 13b for both Magnetic region and geographic region, respectively.This shows the method is likely limited by GOLD's resolution.Recall GOLD pixel size is ∼95 km at nadir, but can stretch to upwards of ∼112 km closer to the limb.Ground based observations see scale sizes on the order of fifty to hundreds of kilometers (e.g., Makela, 2006 and references therein).The separation of EPBs should be considered along with these measures of width as GOLD only sees larger bubbles.GOLD may also see two EPBs that are close together as one.
EPB separations are shown in Figures 12e-12h and 13c-13h.Similar to Section 5.4, further subdividing the data by magnetic region, geographic region, or season shows little change in shape, means, and medians, with the exception of June solstice over SA.Huang et al. (2012) show irregular longitudinal differences between EPBs with data from C/NOFS which appears in Figures 12e-12h and Figures 13c-13h by the wide breadth of separation values possible: ∼200-∼1,200 km for ∼90%-∼92% of the total separations.Quasiperiodic behavior of EPBs observed by Huang et al. (2013) could point to gravity wave seeding with typical separations of 500-1,000 km which is consistent with GOLD's observations.This may explain why the peaks of the distributions, means, and medians consistently lie within the same range, but would require additional gravity wave data to determine.This is not available from GOLD.EPBs were detected and analyzed with coincident GOLD, GNSS TEC, SWARM, ionosonde, and cloud temperature data from 24 October 2018 by Aa et al. (2020) which had interbubble distances of 500-800 km over the SA and AT regions.This shows the existence of a wide range of EPB separations regardless of season or region with a propensity for separations of 500-1,000 km.

Conclusions
An automated method to process GOLD data and detect EPBs has been developed in Section 3 and implemented in Section 4. Simple and robust interpolation methods are used to develop a regular grid (Sections 3.2 and 3.3) that are converted to QD coordinates for analysis (Section 3.4) and combined to form complete north and south images (Section 3.5).Solar and polar influences are removed by implementing a minimum SZA of 107° and QD range of −22° to 22° QDlatitude, respectively (Section 3.6).EPB locations are highlighted via difference images generated by detrending of data with RBAVG (Section 3.7).EIA crest tracking ensures equatorial EPBs are detected (Section 3.8).An algorithm picks and filters out peaks in difference image summations along QDlatitudes corresponding with EPB locations (Section 3.9).These EPB locations are used to calculate a measure of width (Section 3.10) and a final database is created (Section 3.11).At the time of writing, this represents the largest such database of EPBs from airglow observations currently available.
The database is used to select EPBs from distinct images (Section 4.1) to calculate EPB drift speeds (Section 4.2), counts (Section 4.3), and separations (Section 4.5).Double counts are removed and uncertainties in geographic coordinates are computed to ensure unique EPBs from the database are analyzed with reasonable error bounds (Sections 4.4 and 4.6).The results from Section 5 are summarized as follows: • Clear seasonal variation in occurrence of EPBs exists with fewer overall occurrences during the June solstice months and the highest during December solstice.• Slight increases in occurrences can be seen from year to year beyond seasonal variations and are likely due to increased solar activity over the 4 years studied.• EPB occurrences and counts exhibit quasiperiodic behavior in time, regardless of magnetic region.
• Spikes in unique EPBs per day occur during both equinoxes and may surpass December solstice values.
• Overall occurrences by longitude are mostly level with a wide peak west of 35° longitude; however, there is a strong seasonal dependence when split by season.EPB occurrences are more prevalent to the west of GOLD's view during December solstice, are more prevalent to the east of GOLD's view during June solstice, and transition between these two extremes during equinoxes.• EPB drift speeds have a wide range of values, but tend toward ∼100 m s −1 .Further, distributions of drift speeds show little variation by season, geographic region, and magnetic region within GOLD's view.• EPB separations have a wide range of values, but tend toward ∼500 km.Further, distributions of separations show little variation by season, geographic region, and magnetic region within GOLD's view.
When compared to the results from Huang and Roddy (2016), it appears that solar activity has more effect on drift speeds of EPBs than season.Further, C/NOFS data looked at by Huang et al. (2013) show larger EPB separations of 800-1,000 km during solar minimum of 2008 and separations of 400-800 km during 2011.This helps explain the wide range of EPB separations GOLD is able to detect.Taken together and remembering that periods of solar and geomagnetic disturbances have not been removed, season have a strong influence on total detections and unique counts of EPBs within GOLD's purview; however, EPB drift speeds and separations remain largely unaffected by season.This shows the overall behavior of EPBs, within GOLD's purview, is independent of season, geographic region, or magnetic region.Further, winds in the neutral atmosphere are known to change with season (H.Liu et al., 2006).Post sunset ambient zonal ion drifts tend to follow the neutral wind but vary in magnitude (G.Liu et al., 2013).Ion drifts within EPBs can be equivalent to the zonal neutral winds or ambient ion drifts, but they can also differ from either by hundreds of meters per second in the westward direction (Huang et al., 2010;J. Park et al., 2021).This leads to the possibility that the EPBs detected and tracked in this study are not simply drifting at the rate of the background neutral and ion drifts.The nature of how these will look when compared to high solar and geomagnetic activity is the subject for future work.

Figure 1 .
Figure 1.Plots of GOLD data locations from 28 November 2018 at 23:10 UTC in red (left column) and with interpolants overlaid in black (right column).Full view of data (top row), zoomed in (middle row), and nearest neighbor grid used for nearest neighbor interpolation (bottom row).

Figure 2 .
Figure 2. Paired north and south images after (a) Nearest neighbor interpolation and (b) quasi-dipole (QD) conversion for 28 November 2018 at 23:10 UTC.Solar zenith angle and QD cutoffs are displayed and implemented in panels (c, d), respectively.

Figure 3 .
Figure3.Rolling barrel average details and geometry.The pivot point is the current location of the barrel at (x 0 , y 0 ), R 0 is the barrel radius, the region of interest (ROI) consists of all points the barrel could potentially "roll" over, the angles α, β, and θ are used to derive δ for each point within the ROI (refer toPradipta et al. (2015) for further details on these angles).Adapted fromPradipta et al. (2015).
Figure 5a is obtained by implementing the RBAVG along each QD latitude of data from Figure 2d.Finally, RBAVG data from Figure 5a minus data from Figure 2d generates the difference data in Figure 5b.

Figure 4 .
Figure 4. Simple smoothing and rolling barrel average (RBAVG) implemented on a sine function shifted up by 1 (a-b) and the RBAVG implemented on data from Figure 2 along constant QDlatitude of 10 QDdeg (c).Smoothing window is 1.5π with C ≈ 0.707, R 0 in panel (b) is 1.5π, and R 0 for bottom image is four QDdeg.Data is black, detrended data is orange, and difference data is red.

Figure 5 .
Figure 5. (a) Rolling barrel average (RBAVG) implemented on data from Figure 2d, (b) difference data computed by taking RBAVG data and subtracting data from Figure 2d, (c) smoothing Gaussian fit implemented on data from Figure 2d, (d) Gaussian fit line plot along QDlongitude 30 QDdeg of panel (c), (e) quasi-dipole data from Figure 2d with tracked equatorial ionization anomaly (EIA) overlaid in orange, (f) and difference data with tracked EIA overlaid in orange.

Figure 6 .
Figure 6.(a) Quasi-dipole data with summation regions marked in magenta, orange, and red for the northern equatorial ionization anomaly (EIA), magnetic equator, and southern EIA, respectively, (b) raw and filtered sums in red and black, respectively, and (c) uncertainty analysis.Uncertainties about each QDlongitude are ±1.5σ and marked in blue with valid/invalid detections in red/black.

Figure 7 .
Figure 7. (a) Filtered sums with Gaussian fits of peaks overlaid in orange and (b, c) detected Equatorial plasma bubbles marked in red on GOLD disk image pair and quasi-dipole conversion.

Figure 8 .
Figure 8. Image pairs in quasi-dipole coordinates for 28 November 2018 in order: (a-d).Time steps and scan region for pairing are shown above each image as t 0 -N/S, t 1 -N/S pairs (N: north, S: south).Equatorial plasma bubbles (EPBs) are marked in panels (a, b) in red and orange, respectively.Image (a) EPBs are overlaid (c).Plot (e) shows drift speeds as calculated from image (a-c) in meters per second while plot (f) shows EPB separations from image (c).Values are plotted with respect to QDlongitude for quick reference to the EPBs marked above.Error bars are 1 sigma uncertainties derived in Section 4.6.
using image Figure 8a as  0 and image Figure 8c as  1 .Note that a total of 5, 6, and 5 drifts are calculated for the northern, equatorial, and southern regions, respectively, for  1 while there are 7, 6, and 4 EPBs in the same QD regions in  0 .The three new EPBs in  1 (westernmost of the equatorial and southern regions) are assigned the default drift of zero at this stage.EPBs with drift values of zero that are not used to calculate a drift speed are not used in later analysis and are simply used to identify EPBs seen only once.

Figure 9 .
Figure 9.Total unique Equatorial plasma bubble (EPB) detection heat maps by season (a-d) December solstice is from November to February, March equinox is March and April, June solstice is from May to August, and September equinox is September and October.The magnetic equator is not present in heat maps as all EPBs detected in this region are assumed at the magnetic equator.

Figure 10 .
Figure 10.(a-c) Total detections, and (d-f) unique Equatorial plasma bubble counts per night.Solstices and equinoxes have been demarcated by vertical lines and color coded along the top horizontal axis for time series plots: black for December solstice, red for March equinox, blue for June solstice, and pink for September solstice.

Figure 11 .
Figure 11.(a, b) Histograms of Equatorial plasma bubble detections and (c-h) drift speeds by geographic region.Drift speeds are reported in meters per second.Plot (a) shows detections from the entire database whereas (b) shows detections by season.Plots (c-e) and (f-h) show the geographic regions split by magnetic region and season, respectively.Totals and legends are labeled in the top right of each plot with percent of totals along the y-axis.Magnetic regions are "N" for quasi-dipole (QD) northern, "E" for QD equatorial, and "S" for QD southern region.

Figure 12 .
Figure 12. (a-d) Histograms of Equatorial plasma bubble drift speeds and (e-h) separations.Drift speeds are reported in meters per second and separations in kilometers.The full database split by magnetic region are shown in panels (a, e): N for quasi-dipole (QD) northern, E for QD equatorial, and S for QD southern region.The following three rows are the northern, equatorial, and southern magnetic regions split by season.Totals and legends are labeled in the top right of each plot with percent of totals along the y-axis.
Muella et al. (2017) reported average EPB drifts between 96 and 150 m s −1 with a max of 185 m s −1 under the southern EIA over several years during December solstices and equinoxes whileEngland and Immel (2012) used IMAGE-FUV data over the northern EIA to derive drift speeds of 40-160 m s −1 during the December solstices and March equinoxes of 2002 and 2003.These fit well with drift speeds reported in Figures 12b and 12d regarding speeds along the northern EIA and southern EIA by season, respectively.

Figure 13 .
Figure 13.(a, b) Histograms of Equatorial plasma bubble widths and (c-h) separations.Widths and separations are reported in kilometers.Plots (a) and(c-e) display widths and separations over geographic regions, respectively, and subdivided by magnetic region: "N" for northern equatorial ionization anomaly (EIA), "E" for magnetic equator, "S" for southern EIA.Plot (b) shows widths by geographic region.Plots (f-h) subdivide separations over geographic regions by season.
EPB indices start at the western-most region of a given image and move eastward as i increases.Using this notation for the northern regions of Figures8a and 8c would yield Every will therefore have at least one EPB,   , where i is the bubble index within image n and also starts at zero. 1,1 and so on.