Day‐to‐Day Variability of Prereversal Enhancement in the Vertical Ion Drift in Response to Large‐Scale Forcing From the Lower Atmosphere

Equatorial plasma bubbles (EPBs) are plasma density irregularities found in the equatorial ionosphere after sunset and are closely related to the upward E×B drift before it reverses to the generally downward direction at nighttime, referred to as prereversal enhancement (PRE). Although the seasonal variation of PRE and EPB have been generally associated with F and E region dynamo, large variability of EPB occurrence from one day to the next is still puzzling and poses a challenge for equatorial space weather research and forecast. This study shows that the PRE variability is strongly affected by the large‐scale lower atmosphere variability on day‐to‐day weather scales through the E region dynamo at middle latitudes in the summer hemisphere. The whole atmosphere general circulation model with interactive ionosphere electrodynamics employed in the study does not resolve EPB due to insufficient model resolution, but EPB occurrence rate can be deduced from the simulated PRE using an empirical relation. The deduced EPB rate is in good agreement with the observed rates, suggesting the important role of large‐scale dynamics and electrodynamics in preconditioning EPB and the feasibility of probabilistic forecast of EPB using a whole atmosphere model.


Introduction
The rapid ionospheric E region recombination after sunset and the large upward lift of plasma into the F region with reduced recombination rate create a large vertical gradient of plasma density at the equatorial ionospheric bottomside F region at dusk. This large vertical gradient of plasma density is favorable for the development of the gravitational Rayleigh-Taylor instability and thus equatorial plasma bubbles (EPBs) (Kelley, 1989). EPBs grow rapidly as they ascend and elongate along the Earth's magnetic field. The density irregularities associated with EPBs scatter radio signals and can adversely affect communications and Global Positioning System (Kintner et al., 2007(Kintner et al., , 2009. Understanding the cause of EPB, especially the EPB variability, is essential for its forecast. Observations have shown that EPBs are closely correlated with the prereversal enhancement (PRE) of the vertical E × B drift (Anderson et al., 2004;Fejer et al., 1999;Huang, 2018;Huang & Hairston, 2015;Kil et al., 2009). In particular, radar measurements at Jicarmarca radio observatory (JRO) concluded that EPB variability on day-to-day, seasonal and solar cycle scales could be explained mostly by PRE variability (Fejer et al., 1999). Satellite observations further reveal similarities between the longitudinal and seasonal variabilities of PRE and EPB rates (Huang & Hairston, 2015;Kil et al., 2009). However, statistical analysis of the observations also found that the relationship between PRE and EPB is not one-to-one (Huang & Hairston, 2015), and determining other processes affecting EPB occurrence is an active research topic.
Based on the linear theory of the growth rate (Sultan, 1996) and using drift velocity, neutral winds, and ionospheric conductivities from numerical models, longitudinal and seasonal variations of the growth rate have been obtained (Carter et al., 2014;Retterer & Gentile, 2009;Wu, 2015), and they show general agreement with the variations of EPB. However, the day-to-day variation of EPB occurrence, essentially the EPB weather, is still puzzling. Previous studies have suggested various seeding mechanisms for the plasma instability, in particular, by small-scale gravity waves (Fritts et al., 2008;Huang & Kelley, 1996;Hysell et al., 1990;Krall et al., 2013;Kelley et al., 1981;McClure et al., 1998;Retterer & Roddy, 2014). Since gravity waves vary significantly on short time scales, they could be responsible for the daily change of EPBs. However, it was also suggested that the onset of EPB could be spontaneous in the absence of gravity waves (Kudeki et al., 2007).
A distinct seasonal and longitudinal variation pattern has been noted for EPB occurrence rate and PRE. Tsunoda (1985) has found that this seasonal and longitudinal pattern is closely related to the alignment of the geomagnetic field lines and the evening terminator. This is because when the field line and the evening terminator is aligned, the polarization electric charge buildup at the terminator, which is critical for PRE according to the mechanism proposed by Farley et al. (1986), occurs simultaneously over the field line across both hemispheres. As such, the largest integrated effect along the field line is achieved and the PRE is the most intense. Later observational studies have provided clear confirmation of the dependence of both EPB occurrence rate and PRE on the field line-evening terminator alignment (e.g., Gentile et al., 2006;Huang & Hairston, 2015;Kil et al., 2009). Numerical simulations have shown consistent results (Carter et al., 2014;Liu, Bardeen et al., 2018;Retterer & Gentile, 2009;Wu, 2015). These observations and simulations also illustrate, on the other hand, even at time/location of the alignment, the strength of EPB occurrence rate/PRE varies with longitude.
At shorter time scales, the day-to-day variability of PRE from dynamical processes from the lower atmosphere, and its role in EPB variability have received little attention until recently (Ajith et al., 2018;Liu, Bardeen et al., 2018;Carter et al., 2018), as a result of the increasing appreciation that the lower atmospheric forcing can contribute significantly to thermospheric/ionospheric (TI) variability from day-to-day to interannual scales and over broad spatial scales (Liu, 2016). On planetary scales, atmospheric tides are found to affect the plasma structures (Sagawa et al., 2005;Immel et al., 2006), and the semidiurnal tide is in particular important for PRE (Fesen et al., 2000). The variability of these planetary-scale waves, resulting from changes in wave sources and interactions with the variable background atmosphere as well as other atmospheric waves, thus leads to variability of the TI. Such variability are now captured by whole atmosphere models (Akmaev, 2011;Liu, 2014). Of particular relevance to the current study is that the seasonal and longitudinal variability of PRE simulated in the NCAR Whole Atmosphere Community Climate Model with thermosphere/ionosphere extension (WACCM-X) is in good agreement with observations (Liu, Bardeen et al., 2018). The model also reveals large day-to-day variability of PRE.
In this study, the longitudinal, seasonal, and day-to-day variability of PRE are examined using WACCM-X simulations. Building on the current understanding of PRE formation and its seasonal and longitudinal dependence on the field line-evening terminator alignment, we will further examine other factors that may affect the variation of PRE strength over seasonal and longitude. Our analysis confirms the primary importance of F region dynamo in forming PRE and in its seasonal/longitudinal variation and reveals the strong modulation of PRE caused by the E region dynamo of the summer hemisphere. The seasonal and longitudinal variation of F region Pedersen conductivity and PRE strength is influenced by the solar zenith angle and the geomagnetic field strength. The climatological eastward wind in the F region around local dusk leads to the upward E × B drift, while the E region Pedersen and Hall currents associated with climatological winds tend to weaken the upward drift. It is also found that the E region neutral wind is strongly affected by the dynamical variability from the lower atmosphere and plays a critical role in driving the day-to-day variability of the PRE. Although WACCM-X does not have sufficient spatial resolution to directly resolve EPB, daily EPB occurrence rate at each longitude can be deduced using model-simulated PRE and an empirical relation between PRE and EPB occurrence rate (Huang & Hairston, 2015). The monthly average of the deduced EPB rate is shown to agree well with the monthly EPB rate from observations. This agreement suggests the key importance of large-scale dynamics and electrodynamics in driving EPB, as well as the feasibility of probabilistic forecast of EPB using whole atmosphere models.

Numerical Model
WACCM-X is one of the atmosphere components of the NCAR Community Earth System Model (CESM) (Hurrell et al., 2013) (also http://www.cesm.ucar.edu/ for information on the most recent version, CESM version 2). WACCM-X has a top boundary at the upper thermosphere (4.5 × 10 −10 hPa, or ∼600 km). As in the regular configuration of WACCM, the chemistry module is interactive with dynamical transport and exothermic heating (Kinnison et al., 2007). Photochemistry associated with ion species (O + , NO + , O + 2 , N + 2 , N + , and metastable O + states) is part of the chemistry package. Details of an earlier version of WACCM-X can be found in Liu et al. (2010). The recently released version (v.2.0) includes a more self-consistent ionosphere module for WACCM-X that includes computation of electron and ion temperatures, interactive ionospheric electrodynamics including wind dynamo, and O+ transport in the ionospheric F region. In WACCM-X v 2.0, the ion drag and Joule heating are calculated according to Dickinson et al. (1981) and Roble et al. (1982).
The electric field is required by the ion drag and Joule heating calculations. At middle to low latitudes, the electric field is now calculated self-consistently considering forcing by the wind dynamo. The high-latitude electric potential and ion convection patterns are specified according to either Heelis et al. (1982) (parameterized by 3-hourly Kp input) or Weimer (2005) (with 5-min interplanetary magnetic field and solar wind conditions as inputs). The ionization rate, particle precipitation over polar cap and cusp region and neutral heating associated with aurorae are calculated using an analytical auroral model by Roble and Ridley (1987). WACCM-X v 2.0 can be configured either for free running (FR) climate simulations (lower atmosphere unconstrained), or with the lower-middle atmosphere constrained by reanalysis data. The dynamical core of the model has also been recently improved to represent the species dependency of specific heats and mean atmosphere mass in the thermosphere. Validation of WACCM-X v 2.0 shows good agreement with thermospheric and ionospheric observations, including the climatology, short-term variability, and during solar flares and geomagnetically disturbed periods (Liu, Bardeen, et al., 2018;. A FR WACCM-X v 2.0 simulation under solar maximum and geomagnetic quiet conditions is used for the analysis presented in this study. The simulation starts from standard model initial conditions under solar maximum conditions. The solar radio flux at 10.7 cm (F10.7 flux) is set at 200 solar flux unit (sfu), and the geomagnetic Kp index is set at 0.33. The Heelis model (Heelis et al., 1982) is used to specify the high-latitude electric potential and ion convection patterns. Detailed information of the model configuration can be found in Liu, Bardeen et al. (2018). Some of the simulation results, including the composition of neutral and ion species, thermal structures of the thermosphere and ionosphere, seasonal variation of migrating and nonmigrating tides, and ionospheric E × B drifts, are presented in Liu, Bardeen et al. (2018). In this study we will focus on the PRE. The PRE is determined here by finding the maximum vertical E × B drift around the geomagnetic equator between the local times of 17 and 20 hr for each longitude. Vertical E × B drift from only one model level (4.6 × 10 −8 hPa, ∼365 km) is used for the analysis, since it does not change much with altitude in the F region.

Seasonal and Longitudinal Variability of PRE
The equatorial vertical E × B drifts, therefore also PRE, are directly proportional to the polarization electric field in the zonal direction, which is known to be driven by the neutral wind (V) through the ionospheric dynamo process. The electric potential , assumed to be constant along magnetic field lines in the model, is calculated from the dynamo equation (Richmond, 1995) where B is the geomagnetic main field, the electric conductivity tensor, and the conductance is the integration of along magnetic field lines. In general, the winds in the ionospheric F and E regions contribute most significantly to the electric potential, though the E region dynamo is much weaker at night as the E region plasma density decreases rapidly after sunset. However, a detailed analysis of the variability of electric potential in relation to wind and dynamo variability is complicated by global effect of the dynamo terms (the right-hand side) of the potential (elliptic) equation (equation (1)). As shown by Liu and Richmond (2013), wind changes at geomagnetic latitudes up to 50 • can still strongly impact the electric field and drift at equatorial latitudes. In this study, the structure and variability of winds, conductivities and their products (the right-hand side of equation (1)) are examined. Furthermore, time series of equatorial PRE is correlated with time series of the dynamo terms to determine the most impactful processes affecting day-to-day variability. Figure 1 shows the monthly averaged zonal mean zonal and meridional winds (u and v) and Pedersen and Hall conductivities ( P and H ) at local time (LT) 19 hr, when PRE usually reaches its peak. At equinoxes P in the F region ( F P ) is up to 7 × 10 −5 Sm −1 in the Southern Hemisphere (SH) and 4 × 10 −5 Sm −1 in the Northern Hemisphere (NH). This latitudinal structure is similar to results at the same local time from a previous modeling study (Richmond et al., 2015), though the magnitude is larger here due to stronger solar input. F P is larger than the conductivities in the E region ( E P and E H ). At solstices, E P and E H in the summer hemisphere are still large (maxima at 10 −4 Sm −1 and 1.6 × 10 −4 Sm −1 , respectively), and are comparable to or  larger than F P (larger maximum of 10 −4 Sm −1 in the winter hemisphere due to winter anomaly and smaller maximum of 2-4 × 10 −5 Sm −1 in the summer hemisphere), because the summer side E region is still sun-lit or near sun-lit at LT 19 hr in the summer. The zonal wind at LT 19 hr in the F region (above ∼2 × 10 −6 hPa, ∼170 km) is eastward at equinoxes and solstices and becomes westward in the bottomside F region. In the E region the zonal wind at LT 19 hr varies from westward to eastward: In the summer hemisphere the zonal wind is westward at altitude where E P peaks (∼ 2×10 −5 hPa, ∼125 to 130 km) and become eastward at lower altitudes where E H peaks (∼8 × 10 −5 hPa, ∼105 to 110 km). The meridional wind is mostly equatorward in the summer hemisphere at altitude where E P peaks. From Figure 2a, it is seen that the zonal mean zonal wind (average over all local times) in the thermosphere is westward above ∼110 km poleward of 20 • in the summer hemisphere. The mean zonal wind is eastward in the winter hemisphere and in the equatorial region, extending to ∼20 • of the summer hemisphere. The mean zonal wind structure and its hemispheric/seasonal dependence result from the Coriolis torque on the summer-to-winter circulation, similar to the general circulation in the stratosphere. The eastward mean zonal wind around the equator, with wind maximum of more than 20 ms −1 , is often referred to as the superrotation (e.g., Rishbeth, 2002), although the superrotation rate simulated here is weaker than that obtained from observations (Liu et al., 2006). In the upper thermosphere and ionospheric F region, the zonal wind has a strong local time dependence due to the large day-night temperature difference ( Figure 2b) and is always eastward at local dusk in both hemispheres. It is thus opposite to the mean zonal wind in the summer extratropical latitudes at that altitude ( Figure 2a). It is noted that the magnitude and the local time dependence of the zonal wind at the equator shown in Figure 2b is similar to that from the Horizontal Wind Model and somewhat different from the measurements by CHAMP satellite at June solstice (Liu et al., 2006). The directions of the zonal wind at E region altitudes at LT 19 hr (Figures 2c and 2d) are generally the same as the directions of the mean values (thus westward in the summer hemisphere around 125-130 10.1029/2019SW002334 km, and eastward around 105-110 km, also seen in Figure 2a). However, it is seen from Figures 2c-2e that the E region zonal and meridional winds are strongly affected by tidal waves and are thus highly variable, as will be discussed in the next section. The semidiurnal zonal wind perturbations at 1.1 × 10 −5 hPa and 8.4 × 10 −5 are out of phase. This is consistent with the vertical phase of migrating semidiurnal tide, with vertical wavelength between 40 and 50 km.
The formation of PRE results from F region and E region electrodynamics (e.g., Farley et al., 1986;Heelis et al., 1974). According to Farley et al. (1986), the eastward wind in the equatorial F region at local dusk sets up a downward polarization electric field, which maps to the E region along the geomagnetic field lines as an equatorward electric field. This equatorward electric field drives a westward Hall current, and the westward electric current is divergent at the terminator because of the rapid decrease of plasma density (and thus current) on the nightside E region. The current divergence thus builds up negative charges, and along with it eastward (westward) polarization electric field on the west (east) side of the terminator. Building on this understanding, PRE variation could be affected by the following: (i) E region Pedersen current associated with zonal wind (-E P u E B). As can be seen from equation (1) As discussed in section 1, the distinct seasonal and longitudinal variation pattern of PRE is closely related to the alignment of the geomagnetic field lines and the evening terminator. It is also found that even at time/location of the field line/terminator alignment, the strength of EPB occurrence rate/PRE varies with longitude. Such variation is seen in WACCM-X simulations (Liu, Bardeen et al., 2018). Here we explore the factors that can affect this variation. The longitudinal and seasonal dependence of the F region Pedersen conductivity is primarily affected by plasma and neutral densities and temperatures, and geomagnetic field strength (e.g., Heelis, 2004), and the F region plasma density is determined by photoionization before sunset, chemical loss, and transport. If only photochemistry is considered, the ion density is approximately proportional to cos( 2 0 ) ( < 0 ), where is the solar zenith angle, and 0 is the solar zenith angle at the terminator. By further ignoring collision frequency in comparison with ion gyrofrequency, which is proportional to the strength of B, the Pedersen conductivity in the F region is then proportional to cos( 2 0 )∕B 2 (which is referred to factor f S hereafter). is a function of the day of a year, time of the day, and geographic location. 0 depends on the height of interests, and 105.8 • is used (corresponding to 250 km) for our analysis. The f S is calculated for each geographical location at 19 LT, and integrated over 10-50 • geomagnetic latitude in both hemispheres. The combined effect of geomagnetic field and photoionization on longitudinal and seasonal variation can be examined. This latitude range is chosen because the dynamo process between 10 and 50 • geomagnetic latitudes contributes most significantly to the vertical E × B drift at the equator (Liu & Richmond, 2013). This is also seen from the correlation between the PRE and dynamo, as will be discussed later ( Figure 5). The integrated quantity is shown in Figure 3a, and it shows qualitative resemblance to the most prominent longitudinal features around solstices found in observed and simulated PRE (e.g., Huang & Hairston, 2015;Kil et al., 2009;Liu, Bardeen et al., 2018): with largest values located at 40-50 • W in December and January, and two smaller peaks at 0 • and 180 • longitudes around June solstice. These peaks correspond to the locations of the minimum B in the summer hemisphere: the Southern Atlantic Anomaly (SAA) around December solstice and the African and the Pacific sectors around June solstice (Figure 3b). The generally stronger PRE peaks around December solstice than those around June solstice is consistent with the weaker B in the south (especially around SAA). It is also clear that f S is inadequate for equinoctial periods, probably because it does not account for the transport effect, which become increasingly important in the F region, nor chemical loss process.
To facilitate a more detailed analysis of variability on different time scales, zonal and meridional winds, Pedersen and Hall conductivities, and PRE vertical drift are decomposed into the lower frequency (seasonal or longer, denoted by a superscript l) and the high-frequency (denoted by a superscript h) components. This is done by applying a Fourier filter to the time series of these quantities and retaining the components longer (shorter) than 91 days for the high-(low-) frequency components. For multiplication products of conductivities and winds, their low-frequency components are taken as the products of low-frequency components of corresponding terms (i.e., l P U l , l H U l , l P V l ). Their high-frequency components are the total products minus the low-frequency products (e.g., ( P U) h = P U − l P U l ). This is because both the products of the low-frequency components and high-frequency components (e.g., l P U h and h P U l ) and the products of the high frequency components (e.g., h P U h ) contribute to the high-frequency variation of the dynamo process. With this decomposition, the low-frequency variation is analyzed first. The averages of l P U l at 19 LT over  10-50 • geomagnetic latitude in both hemispheres and in the F and E regions are calculated with the vertical ranges of these averages spanning the altitudes where the corresponding conductivities are large (Figure 1). These averages, as shown in Figures 4a and 4b for F and E regions, respectively, are examined as an approximation to the field line integration and thus a measure of the F and E region dynamo forcing (equation (1)). It is evident from the figures that F region dynamo contributes significantly to the eastward electric field (thus upward vertical drift), especially around the two eqinoxes and December solstice, consistent with Heelis et al. (1974) and Farley et al. (1986). The eastward electric field results from the eastward F region zonal wind direction at 19 LT, as discussed earlier, and the seasonal and longitudinal variation is strongly influenced by the variation of P . Although P reaches maximum in December in the summer hemisphere, it is large in both hemispheres during equinox (Figures 1a and 1e), so the integrated effect is large during equinox as seen in Figure 4a. It is noted that this maximum F region dynamo effect at equinox is not seen in Figure 3a. As discussed earlier, f S only considers the plasma dependence on photoionization, but not the transport effect, thus not the two peaks on both sides of the equator. The E region dynamo term ( l P U l ), on the other hand, would on average produce an equatorward current that tends to short-circuit the F to E region current loop, and thus weaken the eastward electric field around both solstices. This is because in the summer side E region where P is still large at 19 hr LT, the zonal wind is mostly westward (Figure 1).
The average E region zonal Hall current and Pedersen current are shown in Figures 4c and 4d. It is noted that the averages of −J E,H x ∕B and −J E,P x ∕B are plotted here, since the eastward zonal current leads to convergence around the evening terminator, and therefore reduces polarization charges and weakens the polarization electric field (thus PRE). In the summer side E region where the Hall conductivity peaks, the mean zonal wind at 19 hr LT is eastward (Figures 1c and 1g). This drives an eastward Hall current, which tends to weaken the PRE. In the E region where the Pedersen conductivity peaks, the mean meridional wind at 19 hr LT is  (Figures 1d and 1h). This again drives an eastward Pedersen current, which likewise tends to weaken the PRE. Therefore, the summer side E region dynamo can significantly modulate the variation of PRE, especially around June and December solstices.

Day-to-Day Variability of PRE and Deduced EPB Rates
As noted in Liu, Bardeen et al. (2018) (their Figure 11 and discussions), the day-to-day variability of PRE is due to model internal variability, most likely variability from meteorological forcing, since the solar and geomagnetic activity conditions do not change in the simulation. Therefore, our analysis will focus on this aspect. Figure 5a shows the correlation coefficients (line contours) between the high-frequency component of P U with the maximum PRE at 155 • E (all at 19 hr LT) around the geomagnetic equator for July. Large positive correlation is noted around 30 • N geomagnetic latitude (∼40 • N geographic latitude) and 1.2 × 10 −5 hPa (∼130 km). Since the standard deviation of ( P U) h is also large therein, the impact of the dynamo is significant. Figures 5b and 5c show the correlation coefficients (line contours) between the high-frequency components of the two E region zonal currents (J E,H x ∕B and J E,P x ∕B displayed in the figures). In regions where the standard deviation of these two currents are large, there are significant negative correlations: ∼ −0.7 for the Hall current at ∼106 km between 20 • N and 30 • N geomagnetic latitude, and ∼ −0.5 for the Pedersen current at ∼115 km around 40 • N geomagnetic latitude. These results clearly demonstrate that the eastward/westward zonal current perturbations indeed tend to weaken/strengthen PRE. The correlation is also evident by comparing the time series of the maximum PRE with those of ( P U) h , −( H U) h , and −( P V) h at locations where the correlations are significant (Figure 5d), in particular during the time period between 17 and 27 July: On 21 July the peak PRE (47 ms −1 ) is much stronger than the climatological value (25 ms −1 ) and corresponds to the maximum of ( P U) h , the maximum of westward Hall current (−J E,H x ∕B), and near the maximum of westward Pedersen current −J E,P x ∕B. All of these quantities drop to their minimum values on 24 July. It should be noted that the zonal wind perturbations at ∼130 and 106 km and meridional wind perturbations at ∼115 km are likely related by a wave (or waves). For example, the zonal wind perturbations at 130 and 106 km would be nearly 180 • out of phase from a migrating semidiurnal tide, and its meridional wind perturbation at 115 km is nearly in phase with the latter. On the other hand, the magnitudes and phases of these waves can be highly variable, so they can cause complex variability of the E region dynamo. It is seen from Figure 5d that the three dynamo terms are in phase on some days, but out of phase on others.
Here we examine in detail the direct cause of large day-to-day variability in the E region over the July period as shown in Figure 5d. We will focus here on the wind variability, and in particular the zonal wind at ∼130 km and 40 • N geographic latitude. Figure 6a is the hourly zonal wind for this time period, with 19 hr LT marked by the vertical lines for each day. While the daily average of the zonal wind is westward and the 19 hr LT zonal wind is often westward, consistent with the summer side E region wind climatology (Figure 1), the daily zonal wind varies significantly from day to day, and at 19 hr LT it can become eastward (e.g., on 21, 22, and 25 July). From the longitude/time plot at this latitude (Figure 6b), it is seen that although the prevailing variation is the migrating semidiurnal tide (denoted as SW2), both the spatial-temporal pattern of the wind and wind magnitude undergo large changes. For example, between 20 and 23 July there is a strong semidiurnal, eastward propagating pattern that contributes to the eastward wind perturbations. This is further quantified by spectral decomposition (longitude and time) of the zonal wind for 19 hr LT (Figure 6c). The spectral reconstruction of zonal wind for zonal wave numbers 0 to 5 and frequency from 0 (daily mean) to 5 cyc d −1 (pentadiurnal) (purple dotted line) is able to capture most of the wind variability. Among them, the semidiurnal eastward components (especially wave numbers 2 and 3, denoted as SE2 and SE3) contribute to the eastward wind perturbations on most days, reaching maximum on 21 July. The semidiurnal westward and daily mean components contribute to the westward perturbations, with the latter being strongest on 24 July. The time series of the reconstructed SE2, SE3, daily zonal mean, and SW2, four of the most prominent components, are shown in Figure 6d. It is seen that the amplitudes of all these components vary from day to day. The phases of SE2 and SE3 also vary significantly, while the SW2 phase is stable. The large variability of the total zonal wind is a combination of the amplitude and phase variability of all these components.
The impact of the E region dynamo on the day-to-day variability is also seen from statistical analysis. Figure 7 shows the monthly standard deviation of peak PRE, as a measure of the magnitude of the day-to-day variability. The standard deviation values are large between June and August at most longitudes except in the American sector. The standard deviation values are also large at multiple locations between November and January. By comparing Figure 7 to Figures 4b-4d, it is seen that the maximum standard deviation of PRE generally corresponds to the time and location with strongest E region dynamo, which is strongly modulated by waves from the lower atmosphere. This is consistent with the case analysis presented in the previous paragraph.
It is seen from Figure 7 that the longitudinal and seasonal variation of the standard deviation is in contrast to the monthly averaged peak PRE values, which are generally weak in the American/Atlantic and Indian sectors between June and August (Huang & Hairston, 2015;Kil et al., 2009;Liu, Bardeen et al., 2018). This suggests the likelihood of having unseaonally large PRE (and thus "abnormal" EPB) in these regions during that time period. This result is consistent with and may provide an explanation to the reported unseasonal development of post-sunset F region irregularities (Ajith et al., 2018;Carter et al., 2018).
From Figure 7 it is also seen that the standard deviation of PRE can also be large around equinoxes (e.g., within the Pacific sector and American sector), when the E region dynamo is generally weak. This may result from F region variability. The ( P U) h in the F region can indeed be large, as seen in Figure 5. The F region dynamo variability is caused by variability of both the zonal wind and conductivity, which is different from the E region. The correlation coefficients from Figure 5 can also be large in the F region, including places where the ( P U) h is large. However, the neutral wind, plasma density (thus electric conductivity), and E × B Figure 8. Monthly EPB occurrence rate, obtained from daily EPB occurrence rate, which is in turn deduced from the maximum PRE using an empirical relation (see text).
drift are closely coupled through transport and ion drag in the F region. It is thus more challenging to disentangle the cause and effect underlying the large correlation between peak PRE and F region dynamo.
From the daily values of peak PRE at each longitude, the probability of bubble occurrence can be deduced using an empirical relation (Huang & Hairston, 2015). Figure 8 shows the monthly averaged EPB occurrence rate deduced from the WACCM-X simulations under solar maximum conditions. In January-February and November-December, the highest occurrence rate (90%) is found in the American/Atlantic sector (60-30 • W). At March and September equinoxes, several peaks are identified: the highest occurrence rate shifts eastward toward the African sector; a smaller peak (∼80%) extends westward into American sector, and another is located in the Pacific sector. Around June solstice, the peak occurrence rate in the African sector drops to ∼60%, and the peak rate in the Pacific sector is ∼70%. The lowest rates are found in the America/Atlantic sector (20%, centered around 70 • W) and in the Indian sector (30%, centered around 70 • E).The deduced occurrence rate is in good agreement with the seasonal/longitudinal variation of the EPB rates obtained from observations (Gentile et al., 2006;Huang, 2017;Huang & Hairston, 2015;Kil et al., 2009). It is noted that the rate values vary among observations, probably due to the different thresholds employed in identifying the plasma depletion (Huang, 2017) and also different altitudes of the measurements. The model-deduced rate values are comparable to those from the ROCSAT-1 (Kil et al., 2009) and C/NOFS (Huang, 2017) measurements, and higher than those from the Defense Meteorological Satellite Program measurements (Gentile et al., 2006).

Discussion and Summary
The ionospheric PRE is shown to vary significantly from one day to the next in WACCM-X simulations. With the solar and geomagnetic forces being held constant, the simulations thus demonstrate that the lower atmosphere is an important driver of the day-to-day variability of PRE. It is known from previous studies that in the mesopause/lower thermosphere and ionospheric E region (MLTE ∼90-150 km), atmospheric waves are large and often reach their peak amplitudes before being damped down by molecular viscosity and ion drag. These waves are also highly variable because of the variability of the wave sources, the background atmosphere they propagate through, and interactions among these waves and between these waves with the large-scale circulation (Liu, 2016). According to the analysis presented above, these waves and their variability can strongly impact the equatorial peak PRE by perturbing the summer side E region dynamo. The E region Pedersen and Hall currents associated with zonal and meridional winds affect the closure of the F and E region current loop and the polarization electric field around the evening terminator. In a specific case study, it is shown that the day-to-day variability of both the amplitudes of phases of migrating and nonmigrating tidal winds are largely responsible for an episode of large PRE variability during July. The importance of summer side E region is also evidenced from statistical analysis: The longitudinal and seasonal variation of the PRE standard deviation is generally similar to the variation of E region dynamo. Therefore, this is likely a key pathway for the lower atmosphere weather to affect the low-latitude space weather. This is further corroborated by the general agreement between the EPB occurrence rate deduced from the simulated E × B drift. While the mesoscale gravity waves may be responsible for directly triggering the EPB, the day-to-day variability of the large-scale dynamics and electrodynamics in the MLTE is probably critical for preconditioning the EPB formation. It should be noted that a recent numerical study by Richmond et al. (2015) found that the low-altitude wind has little effect on PRE. That study examined PRE around equinox, when the E region dynamo is much weaker at dusk. More quantitative comparisons of the dynamo processes at equinox and solstice, as well as further investigation of the day-to-day variability of F region dynamo, are needed in future studies.
This study elucidates the seasonal and longitudinal variation of PRE. The PRE variation is primarily determined by the F region dynamo around dusk and alignment of the field line and the evening terminator (Farley et al., 1986;Heelis et al., 1974;Tsunoda, 1985). While in the F region the monthly averaged zonal wind around dusk is eastward throughout the year, Pedersen conductivity at dusk displays a distinct seasonal and longitudinal variation according to solar zenith angle and geomagnetic field strength. The summer side E region dynamo around dusk can also strongly affect PRE. For example, the climatological westward wind E region wind around dusk tends to weaken PRE by producing an equatorward Pedersen current in the F and E region circuit loop. The peak Hall and Pedersen currents associated with zonal and meridional winds, respectively, are on average eastward and thus tend to decrease the zonal polarization electric field around the evening terminator.
It has been demonstrated in recent studies that, with data assimilation, whole atmosphere models show skills in forecasting the large-scale conditions in the MLTE, including the ionospheric E × B drift (Pedatella et al., 2018;Wang et al., 2014). Therefore, it is possible to use such models to forecast PRE and to determine if the conditions are favorable for EPB formation. Further, by using the statistical relation between PRE and EPB occurrence rate, it is feasible to provide probabilistic forecast of EPB occurrence.
Although EPB occurrence is strongly correlated with PRE, it is not a one-to-one relationship when PRE is between 0 and 40 ms −1 (Huang & Hairston, 2015). This suggests that other processes, large-scale or small-scale, can still be important for EPB. For example, the growth rate of Rayleigh-Taylor instability is also determined by the neutral winds and chemical recombination rates of plasma along the field lines (Sultan, 1996). The temporal change of the vertical E × B drift around dusk may also be important in determining the occurrence of EPB (Chapagain et al., 2009). However, the kind of spatial and temporal coverage needed for exploring these complex processes is not feasible with current thermospheric and ionospheric observations. An approach using a comprehensive whole atmosphere model with interactive dynamics, electrodynamics, and photochemistry that can integrate available observations will be our best hope to investigate and characterize if and how much these factors affect the EPB variability.