© J.M Smith & J. Klenzing, Published by EDP Sciences 2022

1 Introduction

Equatorial plasma irregularities can be observed in both the bottomside and topside F regions of the ionosphere. The prevailing explanation for this behavior is the growth of plasma irregularities in the bottomside into plumes of depleted plasma under the action of a generalized Rayleigh–Taylor (GRT) instability (e.g., Ossakow, 1981; Haerendel & Eccles, 1992; Sultan, 1996). The GRT instability requires some initial perturbation to the bottomside density to initiate the growth into the topside F region. Wave structures in the plasma density produced by traveling ionospheric disturbances (TIDs) and gravity waves propagating up from lower altitudes can serve as seeds. The GRT instability is capable of modeling the plasma plume structures, sometimes referred to as equatorial plasma bubbles (EPBs), that appear similar to those observed by in situ and terrestrial observations (e.g., Keskinen et al., 1998; Retterer & Roddy, 2014; Yokoyama et al., 2014). The prediction of EPBs – including both the occurrence and requirements for a satisfactory prediction – remains an outstanding issue in Space Weather (Carter et al., 2020). Applying the GRT instability calculations based on ionospheric measurements can help to improve predictions of day-to-day variability in EPB occurrence that is not captured in empirical models.

Bubble growth from the bottomside to the topside depends on the strength and duration of the growth rate and the amplitude of initial perturbations. Hypothetically, larger amplitude perturbations can develop into bubbles under the influence of growth rates that are weak, short-lived, or both weak and short-lived (Retterer & Roddy, 2014). Smith & Heelis (2017) found that the spatial scale sizes of the plasma density depletions associated with EPBs, as observed by in situ instruments, were largely constant across several years of disparate solar activity. This suggests that either consistent sources of perturbations for the initiation of the GRT instability or the nonlinear development of EPBs is independent of the initial perturbation. While the spatial sizes are consistent across this time period, the occurrence frequency and depletion depth of bubbles in local time and geographic location can be significantly different for some seasons as well as across the years. Retterer & Roddy (2014) found that the seed perturbations necessary for bubbles are nearly ubiquitous, therefore other large-scale forcing factors likely modulate the bubble occurrence frequency by modifying the ionosphere and thermosphere to enhance or suppress the probability that the seed will form a bubble.

Tsunoda (1985) suggests that the alignment of the magnetic meridian with the dusk terminator is the source of the largest organized seasonal behavior of EPB appearance, as further supported by Burke (2004), Gentile et al. (2006), and Smith & Heelis (2017). On shorter time scales, Carter et al. (2014) suggest that magnetic activity is the primary predictor for EPBs, and Liu et al. (2013) suggest that ultra-fast Kelvin waves propagating from below may also play a role in modulating EPB formation. However, neither of these short-term mechanisms are sufficient to explain the persistent anomalous appearances of bubbles during the June Solstice in the solar minimum period of 2009. Previous work by Smith & Heelis (2018a, 2018b) found that the bubbles observed during this unexpected time were less likely to be observed in groups of discrete plasma depletions extending in local time and longitude than the bubbles observed at typical times, and occurred at later local times although they still had the spatial and density structure of the bubbles typically observed at equinox near dusk. The local time shift in bubble occurrence is accompanied by unique conditions during solar minimum with larger vertical E × B near midnight during 2008/2009 than are typically observed at those local times (Smith & Heelis, 2018a).

Ajith et al. (2016) found that the large growth rates near midnight during the June solstice in 2011–2012 are due to a weak westward electric field in the presence of equatorward meridional wind. This investigation was done for the longitude region surrounding the equatorial atmosphere radar (EAR) at Kototabang, Indonesia at 0.2° S and 100.32° E. In order to investigate the growth rates in this region, C/NOFS drifts averaged over May, June, July, and August of 2011 and 2012 were used to drive the SAMI2 model (Huba et al., 2000). Building on this work, Zhan & Rodrigues (2018) investigated departures in wind patterns from climatological predictions. Zhan & Rodrigues (2018) investigated the equatorial ionosphere near Jicamarca at 11.95° S and 76.87° W during a quiet time event in the June solstice of 2008. Using vertical drifts measured by the Jicamarca observatory to drive SAMI2 they found that weakening downward E × B drifts produced unstable conditions. In addition to neutral winds from HWM93, neutral winds in several configurations were provided to SAMI2 as input. They found that equatorward winds, combined with weakening drifts, can contribute to unstable conditions responsible for post-midnight EPBs. Both Ajith et al. (2016) and Zhan & Rodrigues (2018) are focused on particular longitudes and seasons for comparison and incorporation of ground-based data from Jicamarca and Equatorial Atmosphere Radar (EAR) at Kototabang, Indonesia. Wu (2015) performed a longitudinal and seasonal climatological investigation of the Rayleigh–Taylor instability growth rate using TIEGCM. Investigating the solar moderate year of 2006 with F10.7 ranging from 70 to 100 sfu, they found a seasonal and zonal structure similar to EPB rates observed by DMSP from 1989 to 2004. Their study was primarily restricted to the local time of 1800 h, however positive growth rates near midnight were observed, though not explored and left for future study.

This study presents a seasonal and zonal climatology of GRT growth rates during solar minimum and solar moderate periods under geomagnetically quiet conditions as a function of local time and altitude. The code used to produce these features is released publicly as the growin package, available at https://www.github.com/JonathonMSmith/growin, so that members of the community may, if they choose, run this incarnation of the GRT growth rate calculation for any longitude, season, and year of their choosing for other case studies (Smith & Klenzing, 2020). The growin module was previously used to generate Fejer–Scherliess driven Rayleigh–Taylor growth rates that coincide with bubbles observed by the GOLD mission (Martinis et al., 2021). Here we expand on this by investigating the effect of alternate plasma drift climatologies on the resulting growth rates and compare to the bubble occurrence frequencies observed by the C/NOFS mission as a function of local time.

Other open-source works integral to this project are the pysat and sami2py packages. These are available at https://www.github.com/pysat and https://www.github.com/sami2py respectively. These packages are open-source heliophysics packages that may also be found at https://heliopython.org. This package seeks to adhere to the standards set forth by the heliophysics community outlined by Burrell et al. (2018). Open-source software developed on a free platform like python make the tools for science accessible and transparent to a broader community of scientists. By making the models and data available to readers, comparisons with other models and or data can reveal where future improvements to EPB prediction should be made (Carter et al., 2020).

2 Methodology

The growin packages compute the GRT growth rate from a modeled ionosphere. This is done largely in three steps: growin collates vertical drifts from measurements, provides them to the SAMI2 model, and subsequently computes the flux tube integrated GRT growth rate for the SAMI2 output.

Before SAMI2 is run, the vertical drifts need to be specified. For the purposes of this study, the Fejer–Scherliess drifts (Fejer et al., 1999) and those measured by C/NOFS will be considered. The values for the Fejer–Scherliess model are obtained here from the International Reference Ionosphere (IRI) model. When using drifts from the C/NOFS mission, the drift measurements from C/NOFS must be cleaned and collated. The data from the C/NOFS satellite are collected from CDAWeb using the pysat package (Stoneback et al., 2018, 2019). The drifts (v_⊥m) are in the meridional direction, which is the “upward” direction in the plane of the magnetic meridian, which is perpendicular to B. Drift values are removed where the retarding potential analyzer (RPA) and drift meter quality flags indicate that the plasma density is less than 1000 cm⁻³, drifts are compromised by attitude control, the plasma composition is unacceptable for the RPA to make accurate measurements, or there is possible photoionization on the instrument that compromises the data (Coley et al., 2017). Since we are only interested in geomagnetically quiet periods, drifts are not included when the Kp index has exceeded 3 within the previous 24 h, following Stoneback et al. (2011) and Klenzing et al. (2013). To remove perturbed drifts within the irregularities from the averages, a local average of the plasma density is computed over a 400 km window along the satellite track to exclude points where the normalized deviation from the average is large, ΔN/N > 10%. After the drifts have been cleaned, there are often many gaps in reliable drift data per orbit, therefore drifts are compiled climatologically into seasonal and zonal bins. The zonal bins used here are shown in Figure 1 and collectively include all longitudes. The zones are chosen here to roughly correspond to regions of low variance in magnetic declination. The seasonal bins each consist of three months surrounding the December solstice, March equinox, June solstice, and the September equinox. Even the most sparsely populated season/zone sector contains more than 6000 drift measurements. Within each season/zone sector, a climatological 24 h of local time of drifts are constructed by grouping the drifts further into discreet half-hour bins and computing the median of each bin.

Fig. 1 The yellow lines on this map delineate the zonal boundaries of the four sectors used for this study. The solid red line shows the magnetic equator, while the dashed red line indicates ±10° from the magnetic equator, which is used to restrict vertical plasma drift measurements. These zones are chosen to roughly correspond to regions of low variance in magnetic declination. Drift measurements and bubbles are binned together based on these zonal regions.

The growin module uses the sami2py module to run the SAMI2 model, archive the output, and load the output into python data structures (Klenzing et al., 2019). The sami2py module has been modified from the SAMI2-1.00 version to accept Fourier coefficients as vertical plasma drift input and relies on an updated suite of models. These include NRLMSISe-00 (Picone et al., 2002) to obtain neutral densities and temperature, EUVAC (Richards et al., 1994) to obtain the EUV flux, HWM-14 (Drob et al., 2015) to obtain the neutral wind velocity, and IGRF (Thébault et al., 2015) to obtain the magnetic field geometry and magnitude. For this study, the vertical plasma drifts are specified as Fourier coefficients obtained by fitting the vertical drifts to a Fourier function with ten terms before being given to SAMI2. Both the Fejer–Scherliess drifts from IRI and the C/NOFS drifts are treated in this way. Figure 2 shows two examples of the resulting Fourier fit the C/NOFS data. The black curve shows the median drift values obtained from the cleaned C/NOFS drifts, with error bars showing the average absolute deviation from the median within each bin. The red curve shows the result of the Fourier fit, and the pink dotted line is the linear offset from the fit that is removed when providing the drifts to the SAMI2 model. Finally, the solid blue line indicates 0 m/s. The worst fit to the data is shown on the left and the best fit on the right; these two fits have coefficients of determination (R²) of 0.99 and 0.58, respectively. Table 1 lists the R² values for all of the C/NOFS fits used in this study, including all 16 longitude-season sectors. In order to compare two distinct periods of solar activity during the C/NOFS mission that had high rates of data coverage, the periods of November 2008 through October 2009 and November 2011 through October 2012. These periods are chosen to maintain a contiguous winter solstice. For brevity, these years are often referred to as 2009 and 2012.

Fig. 2 Results of fitting the C/NOFS drifts to Fourier functions. The plot on the left shows the worst fit with an R2 value of 0.58, while the plot on the right shows the best fit with an R2 value of 0.99. The black curve shows the median drift values obtained from the cleaned C/NOFS drifts, with error bars showing the average absolute deviation from the median within each half-hour bin. The red curve shows the result of the Fourier fit, and the pink dotted line is the linear offset from the fit that is removed when providing the drifts to the SAMI2 model. Finally, the solid blue line indicates 0 m/s.

Once the 10 Fourier coefficients are obtained for a particular sector, they can be provided as inputs for SAMI2. SAMI2 is then run for a representative day and longitude for each season and zone for 24 h. Although SAMI2 uses small time steps between 1 and 12 s to compute the plasma motion, for the purposes of this study, we take the data output at approximately 15-min intervals.

Once SAMI2 has completed its run, the output is processed, and flux-tube integrated growth rates are calculated for all flux tubes in each of the 15-min intervals. Equation (1) shows the flux tube integrated linear Rayleigh–Taylor growth rate that is used here. This formulation follows Sultan (1996), where the notation has been updated for consistency:

$$\gamma =\frac{{\sigma }_F}{{\sigma }_{\mathrm{total}}}\left(V-U-\frac{g}{{\nu }_F}\right)K^F-R_T$$(1)

where σ_F and σ_total are the flux tube integrated Pedersen conductivities in the F region, defined here as SAMI2 altitudes above 200 km, and total ionosphere, respectively. V is the flux tube integrated vertical plasma drift or the vertical drift at the magnetic apex. The use of a locally obtained drift value depends on the assumption that the magnetic field lines are equipotentials. U is the flux tube integrated neutral wind perpendicular to the field line weighted by the Pedersen conductivity. ν_F is the flux tube integrated effective F region collision frequency. K^F is the altitudinal gradient of the flux tube integrated electron density at the magnetic equator in the F region. Finally, R_T is the flux tube integrated recombination rate. To calculate all of the flux tube integrated terms, the following quantities, obtained from model and data sources, are required everywhere along an individual magnetic flux tube: magnetic field, plasma drift, electron density, ion density, neutral density, electron temperature, electron and ion gyrofrequencies, electron and ion collision frequencies, and neutral wind.

This approach is built on some common approximations and assumptions. The recombination term from Sultan (1996) is included here, however it should be noted that during nightside local times, the contribution from recombination is negligible due to a minute E region (Huba et al., 1996). Additionally, the formulation for several of the vital quantities are shown in equations (2) through (5).

The Pedersen conductivity is:

$$\sigma =\frac{e}{B}\left[n_e\frac{{\mathrm{\Omega }}_e{\nu }_e}{{\nu }_e^2+{\mathrm{\Omega }}_e^2}+\sum _{i=\mathrm{ionsp}.}\mathrm{}{n}_i\frac{{\mathrm{\Omega }}_i{\nu }_i}{{\nu }_i^2+{\mathrm{\Omega }}_i^2}\right]$$(2)

where n_i/n_e are the ion/electron number densities, Ω_i/Ω_e are the ion/electron gyrofrequencies. The ion–neutral collision frequency, ν_i, is:

$${\nu }_i=2.6\times 10^{-9}\left(n_n+n_i\right)A^{-1/2}.$$(3)

The electron collision frequency, ν_e, is:

$${\nu }_e=5.4\times 10^{-10}{n}_n{T}_e^{1/2}+\left[34+4.18\ln T_e^3/n_e\right]n_e{T}_e^{-3/2},$$(4)

where A is the average neutral density and T_e is the electron temperature. Finally, the vertical gradient in flux tube integrated electron density is:

$$K^F=\frac{1}{N_e}\frac{\mathrm{\partial }{N}_e}{\mathrm{\partial }h}$$(5)

where N_e is the flux tube integrated electron density, and h is height.

Vertical plasma drift as a function of local time at the magnetic equator is obtained from in situ drift measurements made by the C/NOFS CINDI IVM instrument, while the plasma values are obtained from the SAMI2 model, magnetic field magnitude is obtained from IGRF, and neutral wind values are obtained from HWM14. The geomagnetic flux tube is separated into many cells along the flux tube, and all local quantities are calculated for each cell and then integrated. This is done for every field line with an apex altitude in the F region ionosphere. This process is repeated for 24 h of local time at half-hour intervals, again at multiple longitudes representative of the zonal region they occupy, and similarly for each season and year for which vertical plasma drift information is available. This is the finest resolution possible as the vertical drift data used here is only available as a seasonally and longitudinally averaged profile since it is acquired from in situ satellite measurements.

Klenzing et al. (2013) evaluated how SAMI2 responds to the various inputs from C/NOFS and what parameters are best for reproducing quiet time data. This study explored changing the MSIS, EUVAC, neutral wind model, and drift data inputs to the SAMI2 model and focused on the resulting effects on the F-peak and topside ionosphere. They found that the MSIS, EUVAC, and wind model changes modulated the amplitude of the density as a function of local time and the topside density and composition. However, only the use of E × B drifts measured by VEFI onboard the C/NOFS satellite produced the novel decrease of the postsunset F-peak during the solar minimum of 2008. So it would appear that as far as the F-peak and topside ionosphere from SAMI2 is concerned, the most important driver for the changes to the solar minimum ionosphere is driven by the E × B drifts. Figure 3 shows the results of this process using the growin package run using two artificial vertical drift functions. The two panels show the growth rate as a function of local time and altitude driven by a sinusoidal E × B drift function in the left panel and a sinusoidal E × B drift function with a pre-reversal enhancement (PRE) near dusk in the right panel. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The limits are shown in s⁻¹, and the upper and lower limits of ±8e − 4 s⁻¹ correspond to ±2.88 e-folds per hour. An e-fold is the time interval in which an exponentially growing irregularity increases by a factor of e. The vertical drift used as input to SAMI2 is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the altitude corresponding to the maximum value of K^F, or the peak altitude, which will be referred to as the bottomside F region.

Fig. 3 RT growth rates driven by sinusoidal drifts with (right) and without (left) a PRE as a function of local time and altitude. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The vertical drift, V, used as input to SAMI2, is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the peak altitude of KF as a function of local time.

The bubble data set used in Smith & Heelis (2018b) is used here to compute the bubble occurrence frequency as a function of local time. This data set compiles bubbles observed using a rolling ball method to detect density depletions measured by the C/NOFS IVM (Smith & Heelis, 2020). This approach was limited to between 1800 h and 0600 h solar local time, and from 400 km to 600 km altitude, therefore our bubble statistics here are also limited to those regions. To limit our analysis to quiet time bubbles that correspond to the quiet time drifts used to drive the model, bubbles are excluded if the Kp index has exceeded 3 in the previous 24 h, following Stoneback et al. (2011) and Klenzing et al. (2013). The probability of observing a bubble during a particular local time is the number of times that hour contains bubbles divided by the number of times it is sampled. For each satellite orbit, if an hour of local time has more than enough data points to contain at least one bubble of the median width – approximately 200 km – that is considered one sample. Similarly, for each satellite orbit, if an hour of local time contains at least one discreet bubble, it is recorded as such. It is this ratio that is used in the subsequent figures when describing the occurrence of bubbles encountered by the C/NOFS satellite. This process is followed for each zone and season.

3 Results

3.1 C/NOFS during solar minimum

Figure 4 shows the growth rates driven by the C/NOFS climatological drifts in the African sector for all seasons during the 2008/2009 year. The growth rate, vertical drift velocity, and bottomside F region are indicated in Figure 3. The occurrence frequency of bubbles detected by the rolling ball method from Smith & Heelis (2018b) is shown in the purple histograms at the bottom of each panel. During the solar minimum period of 2008/2009, we see that positive C/NOFS driven growth extends up to 2100 h in the September equinox and even occurs near midnight during the June solstice. However, even in the absence of a positive climatological growth rate, there is often still substantial bubble occurrence. This is demonstrated during the March equinox season. In the African March sector, where the geomagnetic geometry favors positive growth, the growth rate is almost exclusively negative from dusk until dawn. Even so, the distribution of bubble occurrence resembles some of those zones and sectors with positive growth rates at dusk and/or near midnight. For example, similar bubble distributions with more strongly positive growth rate profiles appear in the African December sector. The presence of bubbles in the absence of a substantial growth rate may result from coupling to the lower atmosphere, inaccuracies in the model, and physics missing from this calculation of the growth rate. Considering again the alignment of the geomagnetic meridian with the terminator, growth rates in these sectors would be expected to be positive near dusk. Perhaps the growth rates responsible for the observed bubbles are lost in the climatology. However, since the peak bubble occurrence is often in the postmidnight hours, those lost postsunset growth rates resulting from the geomagnetic geometry could only produce bubbles that grow slowly into the topside over several hours. Due to the climatological nature of these analyses, the day-to-day variation of the lower atmosphere from which seeds propagate and the ionospheric environment where bubbles develop cannot be explored here but should be the focus of the future investigation.

Fig. 4 Growth rates driven by C/NOFS drifts function of local time and altitude in the African longitude sector during 2009. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The vertical drift, V, used as input to SAMI2, is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the peak altitude of KF as a function of local time. The occurrence frequency, f of bubbles is shown in the purple histograms at the bottom of each panel.

3.2 C/NOFS during moderate solar activity

Figure 5 shows the growth rates driven by the C/NOFS climatological drifts in the African sector for all seasons during the 2011/2012 year. In all but the June season, we see a pre-reversal enhancement in the drift near dusk, which lifts the bottomside ionosphere. During this uplift, the growth rate is strongly positive along the bottomside. This is consistent with the behavior reported by Tsunoda (1985).

Fig. 5 Growth rates driven by C/NOFS drifts as a function of local time and altitude in the African longitude sector during 2012. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The vertical drift, V, used as input to SAMI2, is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the peak altitude of KF as a function of local time. The occurrence frequency, f, of bubbles is shown in the purple histograms at the bottom of each panel.

There are, however some features of the growth rates driven by C/NOFS drifts that may be surprising. A striking feature seen during the moderate solar conditions during the 2011/2012 year is the presence of large postmidnight and predawn growth rates driven by C/NOFS drifts in the Pacific sector throughout the entire 2011/2012 year coupled with a near absence of a response by the bubble occurrence frequency. The September season is shown in Figure 6 as a representative example of this behavior in the Pacific zone. However, it should be noted that the peak in bubble occurrence frequency typically reaches its maximum at least an hour in local time after the positive growth rate is maximum. These large growth rates being near dawn, are perhaps not given enough time to produce bubbles that rise to the altitude of the C/NOFS satellite. These sectors provide a good location to probe further in a day-to-day study of specific sectors.

Fig. 6 Growth rates driven by C/NOFS drifts as a function of local time and altitude in the Pacific September sector during 2012. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The vertical drift, V, used as input to SAMI2, is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the peak altitude of KF as a function of local time. The occurrence frequency, f, of bubbles is shown in the purple histograms at the bottom of each panel.

3.3 Comparison to Fejer–Scherliess

Figure 7 shows the growth rates driven by Fejer–Scherliess modeled drifts in the African March equinox and June solstice during 2009 and 2012. When inspecting the growth rates produced using the vertical drift from the Fejer–Scherliess model, we can see that in any and all sectors with positive growth rates, that positive growth occurs near dusk in the presence of a PRE. During 2009 in the June sector, the profile of the bubble occurrence frequency suggests positive growth later in the evening. Here the observed CINDI drifts shown in Figure 4 perform markedly better at predicting the presence of near and postmidnight bubbles. However, during June of 2012, the magnitude of growth rates driven by the Fejer–Scherliess drifts are better in agreement with the magnitude of the bubble occurrence frequency than those measured coincidentally with the bubble observations by C/NOFS.

Fig. 7 Growth rates driven by Fejer–Scherliess drifts as a function of local time and altitude during African June and March sectors in 2009 and 2012. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The vertical drift, V, used as input to SAMI2, is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the peak altitude of KF. The occurrence frequency, f, of bubbles is shown in the purple histograms at the bottom of each panel.

There is a bubble occurrence frequency above 0.50 over the African June solstice sector during 2012, when the geomagnetic geometry favors positive growth. This is accompanied by climatological growth rates that are quite small, with a maximum value of 4.16e – 5 s⁻¹. The behavior of the bubble occurrence frequency during the African June solstice sector in 2012 does not appear to be well corroborated with growth rates regardless of the drift inputs used here.

Although the Fejer–Scherliess drift model was developed for measurements at Jicamarca, we have extended the model to all longitudes since it is used as the vertical ion drift component of the IRI model (Fejer et al., 1999; Bilitza, 2018). However, previous work has compared vertical E × B drifts in other longitude zones to the drifts in the Fejer–Scherliess drift model and found multiple differences between the vertical drift behavior of the model and observations (e.g., Oyekola & Kolawole, 2010; Saranya et al., 2014; Yizengaw et al., 2014; Dubazane et al., 2018). So although we expect there to be differences between the Fejer–Scherliess model drifts and the C/NOFS drifts, the key finding is the consistency between the C/NOFS drifts and the bubble occurrence frequency.

As a brief summary of the behavior of the growth rate and its relationship to the behavior of bubble occurrence, Figure 8 provides a broad overview, including all 16 longitude-season sectors. For each sector, the time difference between the occurrence of the maximum growth rate and the maximum bubble occurrence frequency is computed. The distribution of these time differences is displayed in Figure 8. Importantly here, the maximum growth rates are limited to those growth rates occurring between 1700 and 0500 h. The presence of positive growth rates later in local time during 2008/2009 might suggest that the postmidnight bubbles that frequently occurred in 2008/2009 are not the result of slow growth beginning at dusk but more rapid growth beginning at later local times.

Fig. 8 Boxplot showing the distribution of time differences between maximum bubble detection and maximum growth rate. Blue boxes correspond to C/NOFS driven growth rates, while orange boxes correspond to Fejer–Scherliess driven growth rates. The boxes on the left show the distribution during 2008/2009, while the box on the right shows the distribution during 2011/2012. The box extent is one quartile from the median, while the whiskers extend to the fourth quartile. Each black circle marks the time difference for one season/zone sector.

As expected, bubble occurrence frequencies typically increase within a few local time hours of positive growth rates. Furthermore, the growth rates computed using drifts from C/NOFS often perform much better as a predictor for bubble occurrence measured by the same satellite than the growth rates computed from a Fejer–Scherliess driven ionosphere. The distinct solar minimum drifts discovered by many, including Stoneback et al. (2011), produce positive growth preceding the most intense periods of near and post-midnight bubble occurrence. These near and post-midnight bubbles most strongly deviate from the behavior that would be expected as a result of the growth under the influence of Fejer–Sherliess drifts. Sectors with maximum nighttime growth rates below 1e – 4 s⁻¹ typically see the smallest bubble occurrence frequencies during the year.

When considering these results, it is worth bearing in mind at all times the expected season zone dependence on the alignment of the geomagnetic meridian with the terminator. The resulting dependence typically yields favorable conditions for post-sunset bubbles in the African June solstice, African equinoxes, Asian equinoxes, Pacific June solstice, and American December solstice. The effect of this alignment is included in the growth rate calculation via the vertical drift and the total integrated Pedersen conductivity at dusk. This effect is a dusk-time phenomenon since this alignment primarily impacts growth rates via the modulation of the PRE magnitude and does little or none to produce bubbles later in the evening. The expected behavior for bubble occurrence has been observed in the C/NOFS data during moderate solar activity, but during the deep solar minimum of 2008/2009 this expected behavior was not observed Smith & Heelis (2017). The growth rates during 2008/2009, and 2011/2012 also reflect this solar activity dependence on the local time distribution of the bubble occurrence frequency.

4 Discussion

Although the bubble occurrence frequency does appear to react to the growth rate, it does not seem to be true that the strength, or positivity, of the growth rates here, are directly reflected by the bubble occurrence frequency. This could be attributed to seed behavior, as the growth of a bubble depends not only on the strength of the growth rate but also on the amplitude of the seed perturbation (Retterer & Roddy, 2014). The use of the growin package to investigate case studies and day-to-day variability of the growth rate is left for future work.

The growth rates here indicate where a bubble could develop in the presence of a perturbation in density and electric potential. This growth rate defines the rate at which the initial instability develops in the plasma density. It is only after the initial instability forms that plumes can begin to rise into the topside (Sultan, 1996). Once bubbles reach the topside, they may stop rising and remain there until dawn as fossils (Krall et al., 2010). If a fossil bubble has grown high enough to be detected by the spacecraft, we would expect to detect it both during the growth phase and the fossil phase. Since fossilized bubbles have stopped growing, they should not be responsible for the sudden appearance of bubbles later in local time for a given sector. While we do not distinguish between fossilized and growing bubbles in this work, it is noted that the detection of bubbles at later local times may be fossilized versions of the bubbles observed at earlier local times. In this case, we would expect to see bubble occurrence at the same rate later in local time compared to the growth period earlier in the evening.

The growth rates in Figures 4 and 5 are driven by drifts that are measured by the C/NOFS satellite and then averaged in altitude, season, and longitude. Therefore, these climatological results show when the ionosphere is conducive to bubble formation on average. The bubbles discussed here are detected by the C/NOFS satellite between the altitudes of about 400–800 km. Once a bubble forms, the growth rate inside is modified from this background value, and it may continue to grow. This is consistent with the results presented above, showing an early evening region of positive growth followed by a region of negative growth.

Plumes develop and rise at rates that depend on both the amplitude of the seed and the magnitude of the growth rate (Retterer & Roddy, 2014). So under the action of growth rates with the magnitudes seen here, plumes are expected to take some time to develop into the bubbles observed by the C/NOFS satellite. This can be generally seen here as the bubble profiles in sectors with weaker growth rates develop later after a period of maximum growth compared with those regions with large growth rates where the bubble occurrence frequency grows rapidly after the period of maximum growth. Although suppressive growth rates are often seen at local times immediately following conducive growth rates, bubbles continue to appear at those and later local times. Perhaps these are bubbles that formed at lower altitudes during periods of positive growth and are maintained by small positive growth at higher altitudes where the growth rate is not strongly negative.

Sometimes there are negative growth rates at altitudes above strong positive growth, particularly in regions with strong predawn positive growth. Besides its temporal proximity to dawn, this suppressive environment may explain the lack of observed bubbles at the C/NOFS altitudes seen here. Krall et al. (2010) found that when the ion mass density inside the upper edge of a bubble is equal to that of the adjacent background bubbles stop rising. Under this condition, the K^F term of the growth rate would be zero, and as such, the total growth rate would be zero. Therefore under this condition the growth rate calculated, and shown here would be nonexistent, and bubbles stop rising under these conditions. In fact, in those regions where these negative growth rates appear, the K^F term is negative.

The growth rates driven by C/NOFS drifts during 2012 are particularly interesting for their lack of correlation between growth rate and bubble occurrence frequency, an shown in Figure 5. Perhaps the vertical profile of the E × B drift is such that the reduction in E × B drift velocities with altitude is larger than in some other sectors. This is possible since we are not using any model to adjust the E × B with altitude. In the interest of improving growth rate calculations, SAMI2 could be improved upon by modifying the drift model to include the height gradient in the equatorial F region vertical plasma drift (Shidler & Rodrigues, 2019). The limitation here is that these profiles are localized, and a broader collection of measurements from several longitudes would be required to develop a more comprehensive model of the E × B drift in both longitude, season, solar activity, local time, and altitude. For the time being, measured drifts are capable of serving in the place of a global drift model where available.

The Pacific sector, as defined in this study, is much larger than the other zones. Although this zone contains an approximately uniform magnetic declination, the E × B drifts here are still being averaged over a much larger region, and this averaging could produce an overall reduction in the drifts. Perhaps the E × B drifts of the Pacific zone here are more dynamic than the geometry of the geomagnetic field. Additionally, there are daily variations in the winds, the ionosphere, and the lower ionospheric influences that are simply unaccounted for in this climatological approach.

The growin package is intended for long-term use and integration with the modeling and analysis software written, distributed, and maintained by the heliophysics community (e.g., Burrell et al., 2018). Since growin is open source, the methodology implemented in the code is completely transparent and readily modified. Making this code open source and free to download and use means that scientists can explore the code and use it without the need to first get permission or request time and effort on the part of the original developers. This provides an easy avenue for collaboration on the code itself and scientifically by expanding and improving on the capabilities of growin. This lowers the barrier of entry for scientists that want to include the GRT calculation in their study. The work here is a demonstration of some of what can be done when using the SAMI2 model and C/NOFS data to compute the GRT, but the growin package can also be used with any combination of models and data that provide the ionospheric quantities needed to compute the GRT, or even other instabilities utilizing the integration techniques included in this python package.

5 Conclusions

A seasonal and zonal climatology of GRT growth rates during solar minimum and solar moderate conditions as a function of local time was developed here using open source data and software. The Rayleigh–Taylor growth rate calculated using C/NOFS drifts precedes the peak occurrence frequency of bubbles above 400 km altitude, which serves as a better indicator of bubble occurrence than the growth rate calculated using drifts from the Fejer–Scherliess model only during solar minimum. During the solar minimum period of 2008/2009, the C/NOFS drifts produce distributions of growth rates that are far more consistent with bubble occurrence than those growths obtained using Fejer–Scherliess drifts. The growth rates calculated here are consistent with previous studies and affirm that the irregularity occurrence as a function of local time during the deep solar minimum of 2008/2009 is strongly correlated with the unique vertical drift structures under these solar activity conditions. This is particularly true of post-midnight bubbles since Fejer–Scherliess growth rates never indicate positive growth rates after dusk. However, during moderate solar activity Fejer–Scherliess growth rates perform as well as C/NOFS growth rates as indicators of bubble occurrence. And often, the magnitude of the preceding positive growth is a better indication of the magnitude of the bubble occurrence frequency. Certainly, all post-midnight bubbles detected by C/NOFS are not due to later evening growth rates. However, the fact that there are strong climatological growth rates long after dusk may suggest that some, if not many, are produced by conducive conditions later in the evening driven by the distinct positive nighttime drifts of solar minimum.

Provision of the open-source and free-to-use growin module provides easy access to calculation of the GRT growth rate for any equatorial observation. This tool has yielded a broad picture of Rayleigh–Taylor growth rates and provides an accessible platform for further study into the day-to-day behavior of bubble formation and growth. Further study into the growth rate and occurrence of equatorial ionospheric irregularities and bubbles require a hard look at the day-to-day behavior of the ionosphere and the coupled lower atmosphere. A focused exploration of the zones and seasons that have the strongest and weakest agreement between growth rates and bubble occurrence is the best path forward for the growin package developed here. More work is planned to incorporate physical models with increased complexity to investigate other forcing mechanisms, such as the lower atmosphere and magnetosphere.

Acknowledgments

The C/NOFS mission is supported by the Air Force Research Laboratory, the SMC Defense Weather Systems Directorate, the Department of Defense Space Test Program, the National Aeronautics and Space Administration, the Naval Research Laboratory, and The Aerospace Corporation. The solar ionizing flux (F10.7) data are available via the Space Physics Data Facility OMNIWeb interface. The link to this database from OMNI2 is https://omniweb.gsfc.nasa.gov/form/dx1.html. The C/NOFS CINDI data are available from https://cdaweb.sci.gsfc.nasa.gov/ cdaweb/sp_phys/. Work at NASA Goddard Space Flight Center was supported by grant NNH20ZDA001N-LWS and by the Heliophysics Innovation Fund. The study here was generated using pysat version 2.0. Computed quantities utilized output from sami2py version 0.2.0. This work uses the SAMI2 ionosphere model written and developed by the Naval Research Laboratory.The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

Fig. 1 The yellow lines on this map delineate the zonal boundaries of the four sectors used for this study. The solid red line shows the magnetic equator, while the dashed red line indicates ±10° from the magnetic equator, which is used to restrict vertical plasma drift measurements. These zones are chosen to roughly correspond to regions of low variance in magnetic declination. Drift measurements and bubbles are binned together based on these zonal regions.

In the text

Fig. 2 Results of fitting the C/NOFS drifts to Fourier functions. The plot on the left shows the worst fit with an R2 value of 0.58, while the plot on the right shows the best fit with an R2 value of 0.99. The black curve shows the median drift values obtained from the cleaned C/NOFS drifts, with error bars showing the average absolute deviation from the median within each half-hour bin. The red curve shows the result of the Fourier fit, and the pink dotted line is the linear offset from the fit that is removed when providing the drifts to the SAMI2 model. Finally, the solid blue line indicates 0 m/s.

In the text

Fig. 3 RT growth rates driven by sinusoidal drifts with (right) and without (left) a PRE as a function of local time and altitude. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The vertical drift, V, used as input to SAMI2, is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the peak altitude of KF as a function of local time.

In the text

Fig. 4 Growth rates driven by C/NOFS drifts function of local time and altitude in the African longitude sector during 2009. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The vertical drift, V, used as input to SAMI2, is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the peak altitude of KF as a function of local time. The occurrence frequency, f of bubbles is shown in the purple histograms at the bottom of each panel.

In the text

Fig. 5 Growth rates driven by C/NOFS drifts as a function of local time and altitude in the African longitude sector during 2012. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The vertical drift, V, used as input to SAMI2, is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the peak altitude of KF as a function of local time. The occurrence frequency, f, of bubbles is shown in the purple histograms at the bottom of each panel.

In the text

Fig. 6 Growth rates driven by C/NOFS drifts as a function of local time and altitude in the Pacific September sector during 2012. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The vertical drift, V, used as input to SAMI2, is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the peak altitude of KF as a function of local time. The occurrence frequency, f, of bubbles is shown in the purple histograms at the bottom of each panel.

In the text

Fig. 7 Growth rates driven by Fejer–Scherliess drifts as a function of local time and altitude during African June and March sectors in 2009 and 2012. Areas shaded in red and blue correspond to positive and negative growth rates, respectively. The vertical drift, V, used as input to SAMI2, is shown by a solid blue line, and the dashed blue line marks 0 m/s. The dashed black line indicates the peak altitude of KF. The occurrence frequency, f, of bubbles is shown in the purple histograms at the bottom of each panel.

In the text

Fig. 8 Boxplot showing the distribution of time differences between maximum bubble detection and maximum growth rate. Blue boxes correspond to C/NOFS driven growth rates, while orange boxes correspond to Fejer–Scherliess driven growth rates. The boxes on the left show the distribution during 2008/2009, while the box on the right shows the distribution during 2011/2012. The box extent is one quartile from the median, while the whiskers extend to the fourth quartile. Each black circle marks the time difference for one season/zone sector.

In the text