On a simple, data-aided analytic description of the morphology of equatorial F-region zonal plasma drifts

We present results of an effort to evaluate the ability of an analytic model to describe the behavior of the equatorial zonal plasma drifts given inputs provided by readily available climatological models of thermospheric and ionospheric parameters. In a data-model fusion approach, we used vertical drift measurements to drive the model and zonal drift measurements to evaluate its output. Drift measurements were made by the Jicamarca incoherent scatter radar, and model results were evaluated for different seasons and two distinct solar flux conditions. We focused, in particular, on model results for different versions of the Horizontal Wind Model (HWM 97, 07, and 14). We found that, despite its simplicity, the analytic model can reproduce fairly well most of the features in the observed zonal plasma drifts, including the vertical shear associated with the evening plasma vortex. During daytime hours the model predicts similar results for the zonal drifts independently of the HWM used to drive the model. More importantly, the modeled daytime drifts match exceptionally well the behavior and magnitude of the observed drifts for all seasons and solar flux conditions considered. The nighttime results drive the overall performance of the analytic model, and we found that a single HWM cannot provide the best results for all seasons and solar flux conditions. We also examined the main sources of zonal drift variability. Most of the morphology is controlled by the zonal wind dynamo term of the analytic model, but with non-negligible contribution from the vertical drift term. Finally, we examined the contribution from the E- and F-region to the zonal wind dynamo. The morphology of the zonal drifts in the region of observation (240–560 km altitude) is controlled mostly by the F-region winds, but with significant contributions from the daytime E-region particularly during December solstice and low solar flux conditions.


Introduction
Equatorial vertical E × B drifts are well-recognized as playing an important role in the linear stability of the equatorial F-region and in space weather (e.g., Abdu et al. 1983;Fejer et al. 1999;Kudeki and Bhattacharyya 1999;Hysell and Kudeki 2004;Kelley and Retterer 2008;Huang and Hairston 2015;Smith et al. 2015;Huang 2018 role of the vertical drifts in the Generalized Rayleigh-Taylor (GRT) instability (Sultan 1996). The GRT instability is often referred to as being responsible for the wide spectrum of ionospheric irregularities commonly seen in the low-latitude F-region and that are referred to as equatorial spread-F (ESF). Additionally, it has been suggested that the behavior of the zonal plasma drifts as a function of height, and their relative motion with respect to the neutral winds can also play a role in the stability of the equatorial F-region and in the morphology of ESF density perturbations (Kudeki and Bhattacharyya 1999;Hysell and Kudeki 2004;Kudeki et al. 2007;Aveiro and Hysell 2010).
Perhaps more importantly, vertical and zonal plasma drifts provide insight on complex and difficult to directly measure ion-neutral coupling processes taking place at low-latitudes. For instance, previous numerical studies have shown that the behavior of the F-region plasma drifts is heavily controlled by the E-and F-region neutral wind dynamos (Heelis et al. 1974;Eccles 1998;Eccles 2004;Maute et al. 2012;Rodrigues et al. 2012). Observations of thermospheric winds, however, are difficult to make, and current observational techniques are limited to certain times and heights (David et al. 2016). Groundbased Fabry-Pérot Interferometer (FPI) measurements, for example, contribute a significant portion of the observational data used in the development of the empirical Horizontal Wind Model (HWM) (Drob et al. 2015). This technique, however, only provides reliable measurements during clear-sky nighttime conditions and for an altitude of about 250 km (David et al. 2016).
Despite the importance of plasma drifts for the lowlatitude ionosphere-thermosphere (IT) system, groundbased measurements are limited and studies of the shortterm (day-to-day) variability of the drifts are difficult. There is, however, a large number of observations made over a long period (several months to a few tens of years) of time. Long-term observations, such as those made by the Jicamarca ISR, have allowed various climatological studies of the drifts (Fejer et al. 1991;Scherliess and Fejer 1999;Fejer et al. 2005;Hui and Fejer 2015). The limitation in the observations led these studies to focus on heightaveraged values (Fejer and Scherliess 1997;Scherliess and Fejer 1999), despite some case studies having shown significant and important height variations at certain times (Murphy and Heelis 1986;Pingree and Fejer 1987;Kudeki and Bhattacharyya 1999;Fejer et al. 2014) More recently, Jicamarca ISR data sets have become large enough that a successful attempt to estimate the climatological behavior of the drifts as a function of local time and height was possible (Shidler et al. 2019). The observations quantified the average height variations in the drifts, including those associated with the equatorial evening plasma vortex, for different seasons and distinct solar flux conditions. The vortex is the manifestation of a complex ionospheric current system around evening hours, which produces a circular pattern of the ionospheric E × B drifts in the magnetic equatorial plane (Kudeki and Bhattacharyya 1999;Eccles et al. 1999;Rodrigues et al. 2012). A strong height variation (vertical shear) in the zonal plasma drifts is associated with the vortex, with strong eastward drifts being observed at main F-region heights and westward (or weak eastward) drifts in the bottomside F-region. Haerendel et al. (1992) showed that a two-dimensional representation of the low-latitude electrodynamics can describe the height-dependent features of the equatorial electric fields, including the shear in the zonal plasma drifts. In fact, a simplified version of their description is commonly used to explain the observed behavior in the zonal plasma drifts (U i ) and their relationship with thermospheric neutral winds and the vertical component of the plasma drifts (Eccles 1998;Chau and Woodman 2004;Rodrigues et al. 2012;Hui and Fejer 2015;Richmond et al. 2015): where U i represents the magnetic zonal component of the plasma drifts at apex height h. U P φ is the Pedersen conductivity weighted, field-line integrated magnetic zonal wind for a magnetic field line crossing the magnetic equatorial plane at an apex height h. Finally, H and P are the Hall and Pedersen field-line integrated conductivities, respectively, and W i is the vertical plasma drift for apex height h.
Here we present results of a study that investigated the ability of the simple analytic representation (model) of Haerendel et al. (1992) to describe the average behavior of the equatorial zonal plasma drifts given readily available models of thermospheric and ionospheric parameters. We evaluated the contributions from different drivers to the resulting morphology of the zonal drifts. We also quantified the performance of the zonal drift model for different versions of the Horizontal Wind Model (HWM). In a data-model fusion approach, vertical drift measurements made by the Jicamarca ISR are used to drive the model and zonal drift measurements are used to evaluate model performance.
This paper is organized as follows: Section 2 provides information about the analytic model of the equatorial zonal plasma drifts derived by Haerendel et al. (1992). It also provides information about the climatological models and measurements of the thermosphere and ionosphere used as inputs for the analytic model of the drifts. Section 3 presents and discusses the main results of our investigation. Finally, in Section 4, we summarize our main findings.

Analytic model of the zonal plasma drifts
The overall electrodynamic behavior of the equatorial ionosphere is often explained using a two-dimensional (2D) field-line integrated model (Haerendel et al. 1992). In this model, the three-dimensional (3D) current continuity equation is integrated along magnetic field-lines from the base of the ionosphere, where parallel currents are expected to vanish, to produce a 2D current continuity equation in the equatorial plane (Eccles et al. 2015): where R E is the radius of the Earth, L is the radial distance in the equatorial plane measured in Earth radii, and φ is the longitude. The integrated vertical current (J L ), and the integrated zonal current (J φ ) are defined as (Haerendel et al. 1992): where B is the magnitude of the geomagnetic field, E L is the vertical electric field, E φ is the zonal electric field, P and H are the field-line integrated Pedersen conductivity and Hall conductivity, respectively, U P φ and U P L represent the field-line averaged zonal and meridional neutral winds weighted by the Pedersen conductivity, respectively, and U H and U H L represents the field-line averaged zonal and meridional neutral winds weighted by the Hall conductivity. The tilde in˜ P represents the field-line integrated Pedersen conductivity with slightly different geometric factors used in integration. The field-line integrated quantities above are defined in Haerendel et al. (1992) but, for completeness, are reproduced in Appendix B. Additionally, the collision frequencies used in deriving the Pedersen and Hall conductivities are given in Appendix A.
We can rearrange Eq. 3 to obtain an expression for the vertical electric field (E L ): Dividing Eq. 5 by −B and substituting the F-region zonal (U i = −E L /B) and vertical (W i = E φ /B) plasma drifts results in the following expression: In this work, we use Eq. 6 as a simple analytic model for climatological zonal plasma drifts in the equatorial F-region ionosphere. Note that this is simply a more complete version of Eq. 1. The thermospheric and ionospheric parameters used to compute the field-line integrated variables on the right-hand side of Eq. 6 are obtained from widely used empirical models. Geomagnetic field values are determined from the International Geomagnetic Field (IGRF-12) model (Thébault et al. 2015). Neutral densities and temperatures are derived from NRL's Mass Spectrometer and Incoherent Scatter (NRLMSISE-00) model (Picone et al. 2002). Ionospheric densities and temperatures are derived from the International Reference Ionosphere (IRI-2016) model (Bilitza et al. 2017). For this study, neutral winds are determined from three different versions of the Horizontal Wind Model (HWM93, HWM07, and HWM14) (Hedin et al. 1996;Drob et al. 2015). Finally, we determine the climatological variation of the vertical plasma drifts (W i ) from Jicamarca incoherent scatter radar (ISR) observations.
Currently, there are no empirical models specifying values for the integrated vertical current (J L ). However, past studies (e.g. Haerendel et al. 1992;Eccles et al. 2015) examined the potential impact of the last term of Eq. 6 on the morphology of the vertical shear in the zonal drifts during evening hours. They assumed J L to be positive (upward) and roughly constant with height having a magnitude of a few mA/m. The fourth term then predicts westward drifts that are negligible in the topside and stronger in the bottomside F-region where the Pedersen conductance is smaller. This results in an enhancement in the vertical shear of the zonal drifts. While the domain of validity for this assumption still needs to be better investigated, for this analysis we follow previous studies (e.g. Eccles 1998; Chau and Woodman 2004;Rodrigues et al. 2012;Richmond et al. 2015) and neglect the J L term.
All the empirical models used in this study are contained in pyglow, an open source Python module that wraps several upper atmosphere models written in Fortran (https:// github.com/timduly4/pyglow). A brief description of each model is provided in the next few sections, followed by a description of the Jicamarca ISR observations used in this study.

IGRF-12
The IGRF-12 model describes the Earth's large-scale magnetic field beginning in the year 1900 and extends to the present. The magnetic field (B = B(r, θ, φ, t)) is defined in terms of the magnetic scalar potential (B = −∇V ). In IGRF-12, the magnetic scalar potential (V ) is approximated using a finite series expansion in spherical polar coordinates. The best-fit coefficients in the series expansion are determined using measurements of the Earth's magnetic field made by a global collection of groundbased magnetic observatories. A list of the contributing observatories, and a table of the series coefficients for different epochs can be found in Thébault et al. (2015).
In this work, we use IGRF-12 to determine the magnitude of the magnetic field used in evaluating the gyrofrequency used in computing Pedersen and Hall conductivities. We also use the model to trace out the magnetic field lines used in field-line integration. Finally, the spherical vector components of the magnetic field predicted by IGRF-12 are used to determine unit vectors in the magnetic zonal and meridional directions. These unit vectors are used to determine the magnetic zonal (u φ ) and meridional (u p ) neutral winds. A description of the neutral winds used in our modeling is given in the HWM subsection. Shidler and Rodrigues Progress in Earth and Planetary Science (2021) 8:26 Page 4 of 20

NLRMSISE-00
The NRLMSISE-00 is an empirical model of the neutral atmosphere extending from the ground to the exobase, and is an upgraded version of the MSISE-90 model (Picone et al. 2002). It is derived using a combination of ground-, rocket-, and satellite-based measurements to predict neutral temperatures and densities. The model takes as input: location, local time, F 10.7 , and Ap. The model outputs the temperature of the neutral atmosphere (T n ), as well as, number densities for H, He, O, N, Ar, N 2 , and O 2 . In this work, the neutral density and temperatures are used in specifying the ion-neutral collision frequency (ν in ) and electron-neutral collision frequency (ν en ). See Appendix A for additional information.

IRI-2016
The IRI-2016 empirical model is the current version of the International Reference Ionosphere, an international project sponsored by the Committee on Space Research (COSPAR) and the International Union of Radio Science (URSI). The model is developed using measurements of ionospheric parameters made by a worldwide network of ionosodes, incoherent scatter radars, and instruments on-board satellites and rockets.
In addition to location, time, and date, the model takes as input several geophysical indices including F 10.7 , sunspot number, the ionospheric IG index, and the magnetic Ap index. The model predicts several ionospheric parameters including electron density (n e ), electron (T e ), and ion (T i ) temperatures, and ion number densities for H + , O + , He + , NO + , and O + 2 . The height range for IRI-2016 is between 60-2500 km for electron and ion temperatures, and between 75-2000 km for ion composition.
Different model options are available within IRI-2016. In this work, we use the default options except for one. We chose to turn off the storm time prediction for the F2 layer critical frequency (foF2) since we are interested in modeling geomagnetically quiet conditions. Additional information, links to source files, and access to an online version of IRI-2016 can be found in http://www.irimodel. org.
In this work, the ionospheric parameters predicted by IRI-2016 are used to evaluate collision frequencies and the Pedersen and Hall conductivties. See Appendix A for additional information.

HWM93, HWM07, and HWM14
For thermospheric neutral winds we used the Horizonal Wind Model (HWM), which is an empirical model commonly used for global specification of the upper atmospheric general circulation patterns and migrating tides during geomagnetically quiet conditions (Drob et al. 2015). The HWM model makes predictions for the geographic zonal (u, +eastward) and meridional (v, +northward) winds. The model was developed using a large historic data set of the neutral wind measurements inferred from ground-based ISRs, FPIs, and sensors for in situ measurements on rockets and satellites.
Within pyglow, the user has access to three different versions of HWM: HWM93, HWM07, and the most recent HWM14. The primary difference between these versions are the underlying data sets. The data set used to develop the early HWM93 model consists of ∼ 1.2 × 10 6 observations from 35 different instruments spanning a total of 15 years. The data set used to develop HWM07 includes ∼ 60 × 10 6 observations from 35 different instruments spanning over 50 years (Drob et al. 2008). In addition to the data set used to develop HWM07, HWM14 includes ∼ 13 × 10 6 additional observations made between 2004 and 2014 for a total of ∼ 73 × 10 6 observations from 44 different instruments spanning over 60 years (Drob et al. 2015). One major difference between the models is that HWM93 includes the solar flux condition (through the F 10.7 index) as an input, while HWM07 and HWM14 do not. Additional information, including links to the FOR-TRAN source code for HWM, can be found in Drob et al. (2015).
In this work, we use the geographic neutral winds of HWM to determine the magnetic zonal (u φ ) and meridional (u p ) neutral winds used to evaluate the field-line integrated values U P φ , U P L , U H φ , and U H L . See Appendix B for more information.

Climatological plasma drifts from Jicamarca
As mentioned earlier, we take a data-model fusion approach in our analyses. We use measurements of the vertical plasma drifts made by the Jicamarca ISR to aid the implementation of our analytic model of the zonal drifts. Additionally, we use measurements of the zonal plasma drifts made by the Jicamarca ISR to evaluate the accuracy of the model.
The zonal and vertical plasma drifts used in this study were processed as part of a previous investigation (Shidler et al. 2019). In that study, the mean behavior of the zonal and vertical plasma drifts as a function of height and local time were determined from Jicamarca measurements made between 1986 and 2017. The analyses covered an altitude range between 200 km and 600 km for different seasons and different solar flux conditions. The measurements used in deriving the plasma drift profiles are representative of geomagnetically quiet conditions. That is, only observations made when the Kp index at the time and 3 previous values were equal to or less than 3 were considered.
The mean drifts were derived for three different seasons: December solstice (Jan, Feb, Nov, Dec), equinox (Mar, Apr, Sep, Oct), and June solstice (May, Jun, Jul, Aug).

Shidler and Rodrigues Progress in Earth and Planetary Science
(2021) 8:26 Page 5 of 20 Additionally, the mean drifts were determined for two distinct solar flux conditions: low and high. The average F 10.7 of the plasma drift measurements in the low solar flux bin was around 85 SFU. The average F 10.7 in the high solar flux bin was approximately 150 SFU, except during June solstice, where the average F 10.7 was closer to 130 SFU. Of particular importance to this study is that the mean F 10.7 is the same for both zonal and vertical plasma drifts. In their analysis, Shidler et al. (2019) found several cases where the post-sunset plasma drift measurements were contaminated by ESF (Fejer and Kelley 1980). Past studies have found a strong correlation with the occurrence of ESF and the pre-reversal enhancements (PRE) of the drifts (Huang and Hairston 2015;Smith et al. 2015;Huang 2018). The PRE is commonly seen in the evening sector during equinox and December solstice months for the longitude of Jicamarca. Therefore, measurements contaminated by ESF occurred most frequently in the post-sunset sector during equinox, and December solstice, and high solar flux June solstice when the PRE is well developed. Additionally, the data set had few measurements of the zonal plasma drifts in the early morning topside ionosphere when plasma densities are low. During these times, the number of good quality Jicamarca ISR measurements were very limited and mean drifts were not computed. Figure 1 shows the JRO vertical drifts used as input to the zonal drift model. The top (bottom) row shows the vertical drifts for each season during LSF (HSF) conditions. Regions in black indicate times and heights where drift measurements were limited in number and did not produce reliable statistics. As mentioned above, most of the limitation in measurements were caused by ESF coherent echoes contaminating the ISR measurements and by low plasma densities causing ISR echoes with low SNR. The JRO zonal drifts will be shown along with modeled zonal drifts in the following sections.

Results and discussion
Figures 2 and 3 summarize the behavior of the climatological zonal drifts estimated from Jicamarca observations and the zonal drifts predicted by the analytic model for low solar flux (LSF) and high solar flux (HSF) conditions, respectively.
Results for each season are grouped in the columns for each figure. The top row shows the climatological zonal drifts observed by the Jicamarca ISR. The second, third, and fourth rows show the results for the modeled zonal drifts using HWM14, HWM07, and HWM93 as neutral wind input drivers, respectively. For the LSF conditions, the geophysical indices used by the driving empirical models were F 10.7 = 85 SFU and Kp = 2. For HSF conditions, we used Kp = 2 and F 10.7 = 150 SFU during Equinox and December solstice, and F 10.7 = 130 SFU during June solstice. Day-of-year (DOY) numbers 80, 172, and 355 were used to represent equinox, June solstice, and December solstice, respectively.

On the behavior of the observed zonal drifts
The observed zonal drifts (top panels in Figs. 2 and 3) were discussed in detail by Shidler et al. (2019). Here, we limit ourselves to summarize some of the  The results also show that height variations in the zonal drifts are negligible during daytime hours. At night, however, height gradients in the zonal drifts are noticeable, particularly in the evening sector, between about 1800 LT and 2000 LT at bottomside F-region heights (below ∼350 km). This height variation is associated with the vertical shear in the zonal plasma drifts that is part of the evening plasma vortex (Kudeki and Bhattacharyya 1999;Rodrigues et al. 2012). The drift measurements are not as severely affected by ESF during June Solstice and the shear can be better identified in the observations for the period. Nevertheless, despite data gaps, the evening height gradients in the zonal drifts can also be observed in December solstice and Equinox. Finally, the observations in Figs. 2 and 3 also show that the evening shear in the zonal drifts is more pronounced during HSF conditions.

On the behavior of modeled zonal drifts
We now turn our attention to the modeled zonal plasma drifts shown in the lower three rows of panels in Figs. 2 and 3. Again, results were obtained for three different HWM models as inputs to our analytic model.
The most striking result is that the simple analytic model is capable of producing most of the features seen in the observed zonal drifts. As expected, the analytic model shows the predominant diurnal variability of the drifts, being westward during the day and eastward at night. The model results also show that the amplitude of the zonal drifts is larger during nighttime compared to daytime values. Additionally, in agreement with the observations, the modeled drifts show negligible variation with height during daytime hours but significant height gradients at night, especially in the evening sector during HSF conditions. Nevertheless, while the overall features of the observed drifts are reproduced by the analytic model, more detailed analyses show quantitative differences in the model results with respect to the observations. These differences are more explicit during nighttime hours when the drifts show more variability with height. The magnitude of these differences vary with the HWM model used.
In order to aid our analysis and comparison of the drifts, Figs. 4 and 5 show height profiles of the measured and modeled drifts for LSF and HSF conditions, respectively. The profiles also serve to show the variability of the observed drifts used in the averages, and help us to assess how the models compare to the measurements. The black solid lines represent the observed (average) drift profiles. The horizontal bars represent the variability (standard deviation) of the measurements used to compute the averages. Solid lines of red, green, and blue represent model results using HWM93, HWM07, and HWM14, respectively.
The profiles show that the variability of the observed zonal drifts is larger during nighttime hours. This is, at least in part, due to larger errors in the measured drifts associated with ISR observations made under reduced nighttime electron densities (lower SNR echoes). Considering all the observations, the median variabilities for daytime (0600 -1800 LT) and nighttime (1800-0600 LT) hours were 29 m/s and 56 m/s. In some extreme cases, however, particularly in the topside during nighttime LSF conditions the variability exceeded 90 m/s. More importantly, the results in Figs. 4 and 5 show that the analytic model tends to predict zonal drifts within the variability of the observations during most of the local times. This is true for all the three HWM models.
Another noteworthy finding that can be better observed in the drift profiles is the convergence of all model results to very similar drift values during daytime hours, from about 0800 LT to 1600 LT, for all seasons and solar flux conditions. The convergence of our analytic model for all the HWM drivers during daytime hours is a result of the relatively weak differences in the field-line averaged zonal neutral winds (U P φ ) between each HWM model. During the daytime, differences between U P φ for each of the different HWM wind models are typically around 5 m/s to 20 m/s, while at night the differences can exceed 80 m/s. Finally, the results in Figs. 4 and 5 show a remarkable agreement between modeled and observed drifts during daytime, especially considering the reduced error and variability of the drifts during that period. The accuracy of physics-based models of the equatorial plasma drifts, including our simple model of the zonal drifts (Eq. 6), rely heavily on predictions of the neutral winds. This is because the neutral wind dynamo is the primary driver for the polarization electric fields and associated E × B plasma drifts observed in the equatorial F-region ionosphere (Heelis et al. 1974;Eccles 1998;Eccles 2004;Heelis 2004;Maute et al. 2012). Consequently, our results serve as an independent indicator of a good representation of the daytime neutral winds by the HWM models and of the daytime conductivities.
An overview of the instrumentation and observations used to derive daytime predictions of the meridional and zonal winds in HWM93 are given by Hedin et al. (1988), Hedin et al. (1991), and Hedin et al. (1996), but sparse spatial and temporal coverage of the data resulted in several technical challenges in constructing this model (Drob et al. 2008). Daytime model predictions of the neutral winds in the more recent HWM07 and HWM14 versions benefited from a large increase (∼ 55 × 10 6 data points) in neutral wind observations made by two instruments on the Upper Atmospheric Research Satellite (UARS) (Drob et al. 2008) The first instrument is the Wind Imaging Interferometer (WINDII) which uses Doppler shifts coming from 557.7 nm O( 1 S) optical emission lines to measure Fig. 4 Altitude profiles for the three seasons during LSF conditions. The blue, green, and red lines represent model zonal plasma drifts using HWM14, HWM07, and HWM93, respectively. The black curves and error bars represent the average and standard deviations of the Jicamarca zonal plasma drifts, respectively   (Hays et al. 1993). Therefore, the results in Figs. 4 and 5 serve to confirm the reliability of the neutral wind observations made by the instruments mentioned above to predict realistic daytime zonal drifts at Jicamarca. We quantify this assessment in the Section 3.5 by computing the RMSE for the model predictions with respect to the observations. Finally, we must comment on the ability of the analytic model to reproduce the vertical shears in the zonal drifts, which are observed during evening hours. Looking at the profiles for the evening sector (1800-2400 LT) in Figs. 4 and 5, one can see that model predictions do not perfectly match the observed average values but are well within the variability margins in most cases. The results also show similar drift predictions for models driven by HWM07 and HWM14 especially during LSF conditions, which indicates similar wind predictions for these models. A noticeable departure from good model-observation agreement, however, occurs for HMW93 profiles between 1900 LT and 2200 LT during LSF conditions (all seasons), and between 1800 LT and 2100 LT for HSF conditions (June solstice). The HWM93-driven drifts tend to overestimate the observations. This departure is most likely caused by HWM93 overestimating the eastward zonal neutral winds, and by extension the first term in Eq. 6, in the F-region during these times. Monthly comparisons between the three HWM models at 250 km (at a constant F 10.7 value of 107 SFU) with FPI and satellite observations of the zonal winds found that HWM93 overestimates the data by more that 50 m/s just after dusk between 1800 and 2100 LT (Drob et al. 2015). A second obvious departure from the overall model-observation agreement occurs for HWM14 and HWM07 profiles between 2100 LT and 2400 LT during HSF conditions (December solstice and equinox). The drifts derived using these models tend to underestimate the observed drifts. A possible explanation for these results is that HWM14 and HWM07 do not have a solar flux dependence in the predictions. Additional details will be provided in the Section 3.5.

On the sources of the modeled drift morphology
The previous section showed that the analytic model is capable of reproducing most of the features of the zonal drifts independent of the HWM model used to drive it. The magnitude of the modeled drifts are, in most cases, within the variability of the observations. Disagreements between models and observations can be found, however, and occur mostly in the evening sector, a period of significant changes in conductivity and in the relative contributions from E-region and F-region dynamos.
In order to better understand the drivers of zonal drift morphology, we now take a closer look at the contribution of the different terms in Eq. 6 and how they control the behavior of the zonal drifts. For this analysis, we start by looking at June solstice HSF conditions when data is available for most local times and heights.
The top row of panels in Fig. 6 show the resulting modeled zonal drifts for each HWM model input. The left-hand side panels of Fig. 6 show results for HWM93, the middle panels show results for HWM07, and the righthand side panels show results for HWM14. The other row of panels show the contributions from the three first terms in the right-hand side of Eq. 6. The terms are indicated on top of each panel. Note the different scales on the colorbar for each term. We remind the reader that the analytic model neglects the contribution of the J L /(B P ) term.
The second row (from the top) of panels in Fig. 6 show that most of the morphology of the analytic model zonal drifts (U i ), particularly the diurnal variation and the post-sunset vertical shear, is heavily controlled by the first term in Eq. 6, that is, by the winds in the magnetic zonal direction weighted by Pedersen conductivity (U P φ ). This behavior was also observed by Rodrigues et al. (2012), who performed similar term-comparison analysis between 1700 LT to 2200 LT using TIE-GCM model run. Eccles (1998) also provides reasons for neglecting all but the first term (U P φ ) of Eq. 6 in the development of their simple model of the zonal plasma drifts.
The third row of panels show that the second term of Eq. 6, associated with magnetic meridional winds, has negligible contribution to the morphology of the zonal drifts. In the original formulation of the simple model of the zonal drifts, Haerendel et al. (1992) neglected contributions coming from meridional winds stating that they are small at near-equator latitudes during the evening. Other studies also followed this formulation, and neglected meridional wind contributions (Rodrigues et al. 2012;Richmond et al. 2015). Eccles (1998) used the 2D field-line integrated approach to develop a simple model of the zonal drifts to study the behavior of the low-latitude electric field throughout the day and at different altitudes. They make the claim that U H L and W i have the same sign and magnitude throughout the day, and therefore, the second and third term of Eq. 6 cancel each other.
For completeness and verification, we included the contribution from the Hall-weighted meridional winds (U H L ). The third row of panels in Fig. 6 show that U H L and W i , in general, have the same sign, but U H L has a magnitude of only a few m/s throughout the day. We find similar behavior for all seasons and solar flux conditions. We must point out that while meridional winds predicted by models do not seem to contribute significantly to the morphology of the zonal drifts, additional efforts have yet to be devoted to a better understanding of the accuracy of the wind model predictions. For instance, the recent study of Navarro and Fejer (2019) found a general agreement between quiet-time averages of the nighttime northward wind observations from three FPI instrument in Peru with HWM14 model predictions during moderate solar flux conditions. They noted, however, that HWM14 underestimated the early night southward winds during equinox, and significantly overestimated northward winds in the post-midnight sector during December solstice.

Shidler and Rodrigues Progress in Earth and Planetary Science
Finally, the bottom row of panels show the contributions from the third term of Eq. 6, which is independent of the wind models. This term is controlled by the product between the vertical drifts and the Hall-to-Pedersen conductivity ratio. The results show that the contribution from the third term is non-negligible during daytime hours (up to 15 m/s) at all F-region heights. At night, the contribution is only noticeable (up to 5 m/s) during evening hours (1900 -2300 LT) at bottomside F-region heights.

On the contribution of H P W i : daytime versus nighttime
The contribution of the third term in Eq. 6 to the morphology of the zonal drifts during daytime hours comes from enhanced Hall conductivities and increased Hallto-Pedersen conductivity ratios. In the evening, the ratio is only non-negligible for a few hours after sunset and for shorter magnetic field lines, that is, for field lines with apex heights below the equatorial F-region peak. For instance, during June solstice HSF conditions, the ratio of the field-line integrated Hall-to-Pedersen conductivity ( H / P ) is ∼0.8 during the day and ∼0.2 at night, except at bottomside F-region heights during evening hours when H / P ∼ 0.6 at an apex height of 240 km. Therefore, the contribution to the zonal drifts by the third term in Eq. 6 during the day is on the order of the daytime upward drifts. As a result, the third term contributes significantly to the behavior of the daytime zonal drifts. At night, the contribution is only a fraction of the magnitude of vertical drifts, except at lower apex heights. As a result, the evening downward drifts contribute to eastward drifts in the bottomside F-region weakening the shear in the zonal drifts. Figure 7 now shows the results for June Solstice LSF conditions. The results are similar to those seen for HSF conditions, that is, the first term on the right side of Eq. 6 controls most of the morphology of the zonal drifts. Again, the contribution from the second term is negligible. The contribution of the third term is even more noticeable during LSF than what it was for HSF. It reaches up to 20 m/s during the day and 15 m/s at bottomside F-region heights during nighttime hours.  Figure 1 shows that the nighttime downward drifts are weaker during June LSF compared to June HSF, and so the enhancement in the bottomside eastward drifts predicted by the third term of Eq. 6 is driven by a strong increase in the H / P ratio for the same region of apex heights. Figure 8 shows the local time dependence of the Pedersen conductance (left), Hall conductance (middle), and Hall-to-Pedersen ratio (right) for June solstice and at an apex height of 240 km. The solid and dashed lines correspond to HSF and LSF conditions, respectively. Figure 8 shows that the Hall conductance is approximately constant throughtout the night, and decreases from ∼4.0 mho during June HSF to ∼2.6 mho during June LSF. Therefore, the local time behavior of the Hall-to-Pedersen ratio is controlled by the Pedersen conductance. At 2000 LT, the Pedersen conductance for both LSF and HSF is ∼7.5 mho and the Hall-to-Pedersen ratio reaches its maximum nighttime value during HSF. However, after 2000 LT the Pedersen conductivity increases (decreases) during HSF (LSF), and the corresponding Hall-to-Pedersen ratio decreases (increases) until the early morning. The different behavior of the Pedersen conductance during HSF compared to LSF conditions is driven mainly by a decrease in the nighttime F-region plasma densities occurring during LSF conditions.

On the contribution of H P W i during equinox
We now turn our attention to Equinox conditions when the vertical drifts are known to maximize during evening hours as a result of the PRE (Fejer et al. 1991). The results for HSF and LSF conditions are shown in Figs. 9 and 10, respectively. The presentation of equinox results follows the same format of Figs. 6 and 7 for June solstice conditions. Overall, they show similar relative contributions from the three terms of Eq. 6 just discussed for June solstice. However, a couple of interesting differences can be noticed in the third term and are discussed here. The results for equinox HSF conditions (Fig. 9), show, despite reduced data due to ESF, that the third term in Eq. 6 (bottom panels associated with W i ) contributes even more significantly to the shear in the zonal plasma drifts, especially between 1600 and 2000 LT, than previously found for June solstice. This is a result of the large upward drifts and enhanced H / P at bottomside Fregion heights in the evening sector seen during equinox in the Jicamarca sector.
In addition to a direct effect on the third term of Eq. 6, the enhanced vertical drifts also have an indirect effect on the zonal drift model. During periods of large PRE, the F-region is expected to be lifted up to higher altitudes. This is what is seen in the IRI-2016 predictions, and it serves to show a good agreement between the empirical models of ionospheric density and drifts. As large vertical drifts move the F-region upward, it lowers the P at The results in Fig. 10 show that this effect is not as pronounced during Equinox LSF conditions. This is a result of the severe reduction in the PRE peak (Fejer et al. 1991) during conditions of low solar activity. Another interesting feature observed in the zonal drifts during equinox is the enhanced contribution of the vertical drift (third) term to the modeled zonal drifts in the pre-sunrise sector during LSF conditions (Fig. 10).
The bottom panels of Fig. 10 show an enhancement in the eastward drifts around 0330 LT, particularly at lower apex heights when eastward drifts can be as large as ∼20 m/s. Similar to our previous analysis (see Fig. 8), an enhancement can be seen in the H / P ratio during equinox around this time, but that ratio is ∼1.25 compared with ∼1.75 during June solstice. This is mainly due to the empirical models predicting larger pre-sunrise Fregion Pedersen conductivities during equinox compared with June solstice. However, Fig. 1 shows that the presunrise downward vertical drifts are larger during equinox

Contributions from E and F regions to U P φ
We already pointed out that a significant portion of the zonal drift morphology comes from the zonal wind dynamo term (U P φ ) of Eq. 6. Here, we examine the contributions from the E and F regions to the dynamo as predicted by the empirical models for the Jicamarca location. Figures 11 and 12 show the results of our estimates for LSF and HSF conditions, respectively. The top panels in each figure show the total zonal wind dynamo (U P φ ) contribution to the zonal drift model. For this analysis, we only used results driven by the HWM14 model. Each column represents a season. The middle panels represent the contribution from the F-region (U PF φ ), that is, from points along the magnetic field lines that are above 150 km altitude. The bottom panels, on the other hand, represent the contribution from the E-region (U PE φ ), that is, from points along the magnetic field lines that are below 150 km altitude. Note the different scales in the colorbars.
The results show that the F-region controls most of the zonal wind dynamo term during nighttime LSF and HSF conditions for equatorial altitudes above 240 km. This is in good agreement with previous observational results that indicated strong neutral-plasma coupling at F-region heights during nighttime. For instance, several studies have shown good agreement between zonal plasma drifts measured by ISRs and zonal neutral winds measured by FPIs (Biondi et al. 1988;Fejer 1993;Chapagain et al. 2013;Navarro and Fejer 2019;Navarro and Fejer 2020) or by sensors on satellites (Fejer et al. 1985). Furthermore, our results are in good agreement with studies that showed that low-latitude plasma drifts and their coupling to neutral winds is controlled not only by local parameters but by the distribution of winds and conductivities along magnetic field lines (Coley et al. 1994).
The E-region contribution to equatorial F-region zonal drifts above 240 km altitude is noticeable and must be considered during daytime hours. This is particularly true for LSF conditions as shown in Fig. 11. We found that the E-region dynamo contribution to daytime F-region zonal drifts decreases with solar flux as indicated by the results presented in the bottom panels of Figs. 11 (LSF) and 12 (HSF). Our results also show that the absolute E-region contribution to daytime zonal drifts is minimum during June solstice, increases during equinox and maximizes during December solstice.
We also found that our model only shows a significant contribution of the E-region to the nighttime zonal drifts during December solstice between about 0000 and 0600 LT. This contribution has magnitudes reaching up to 20 m/s near sunrise during LSF. In fact, the E-region contribution is comparable to the F-region contribution at that time. Our analysis showed that during this time, HWM14 predicts strong westward winds (up to 65 m/s) at points near the E-region conductivity peaks, especially south of the magnetic equator. These points make a significant (2021) 8:26 Page 15 of 20 Fig. 11 The contributions from the E and F regions to the zonal wind dynamo during low solar flux (LSF) conditions. HWM14 was used to drive the zonal drift model contribution to the westward zonal drifts predicted by the E-region dynamo field-line integral (U PE φ ). For those same field line points, HWM14 predicts weaker westward or even eastward winds during equinox and June solstice.
Finally, our results show that the empirical models do not predict any significant contribution from the E-region to the vertical shear in the evening zonal drifts for the altitude range covered by our observations, that is, between 240 and 560 km at the magnetic equator. Using results of a self-consistent numerical model (TIE-GCM), however, Rodrigues et al. (2012) found that the E-region can contribute to the shear. Their analyses, however, extended to altitudes as low as 200 km and their results show significant E-region contributions below about 250 km at 2000 LT during solar flux conditions representing an F 10.7 = 220 SFU.

On model performance
In this section we take a more quantitative look at how the analytic model, driven by different versions of HWM, compares to observations.
In order to make a quantitative assessment we used two metrics, the root-mean-square error -RMSE (σ = i w i (y i −ŷ i ) 2 / i w i ) and the model bias (μ = i w i (y i −ŷ i )/ i w i ). Here, y i represents an observed zonal drift value,ŷ i , represents a modeled zonal drift value, and w i is the weight and is given by the inverse of the standard deviation for measurements used to derive the average zonal drifts. These metrics were computed for each season and solar flux condition. Table 1 shows our results of the RMSE and model biases for the full day (0000-2400 LT), daytime hours (0600 -1800 LT), and nighttime hours (1800-0600 LT). The best performing model, using the RMSE as a metric, is underlined. The most important result in Table 1 is that a single HWM model cannot produce the best results (lowest RMSE values) for all the seasons and solar flux conditions. Table 1 also shows that all of the models perform with reduced errors during daytime compared to nighttime. There are a few cases, such as HWM07 results for June Solstice and HWM93 results for LSF December Solstice, which show similar performance for nighttime and daytime conditions. This quantitative assessment of the model performance confirms our inference about the accuracy of the daytime neutral winds predicted by the HWM models as well as the empirical model estimates of the ionospheric conductivities.
Additionally, Table 1 shows that the performance of the model driven by different HWM varies with solar flux conditions. At LSF conditions, the best performances (full day RMSE) are shared between model runs driven by HWM07 and HWM14. At HSF conditions, however, the model runs driven by HWM93 outperform the other two model runs, except for June solstice where HWM14 is the best performing model.  These results can be explained, to a large extent, by the fact that the majority of the observations used to develop HWM14 were made during periods of moderate-to-low solar activity (average F 10.7 of 107 SFU) (Drob et al. 2015). The better performance of HWM93 during December solstice and equinox HSF conditions can be attributed to a reduced number of HSF observations used to develop HWM07 and HMW14 and to the solar flux parameterization used by Drob et al. (2008) and Drob et al. (2015). The best peformance by HMW14 during June solstice HSF can be explained by the reduced number of Jicamarca ISR observations for high solar activity during June solstice (see the Section 2.6). The observations for June HSF had a mean F 10.7 of 130 SFU instead of 150 SFU for equinox and December solstice. Therefore, the June HSF observations used to evaluate the model better reflect moderate solar flux conditions, which are expected to be better represented by HWM14 than HWM93.
We now take a look at the model biases (μ) provided in Table 1. This metric serves to evaluate whether the model zonal drifts overestimate/underestimate the observed average zonal drifts. Focusing on nighttime, when model predictions and average zonal drifts are eastward (positive), we find that for most seasons and solar flux conditions, HWM93 tends to overestimate the average zonal drifts, while both HWM14 and HWM07 tend to underestimate the average zonal drifts. Given the importance of the nighttime F-region wind dynamo to the morphology of the drifts it is worth mentioning that Navarro and Fejer (2019) also found that HWM14 tends to underestimate the quiet-time zonal winds in the pre-midnight sector during moderate solar flux conditions. During daytime, when model predictions and average zonal drifts are westward (negative), all three models tend to overestimate the average zonal drifts. However, for a given season and solar flux condition, model biases during daytime hours have comparable values further indicating that all three HWM models tend to converge during day.

Conclusions
Currently, a climatological model of the equatorial ionospheric zonal E×B drifts as a function of both local time and height has yet to be developed. While height variations are known to occur especially in the evening sector (Kudeki and Bhattacharyya 1999;Rodrigues et al. 2012;Richmond et al. 2015;Shidler et al. 2019), the lack of adequate ground-based radar observations forced previous studies to focus on height-averaged zonal drifts (Fejer et al. 1991;Fejer et al. 2005).
Here, we present results of a new effort that examined the ability of widely available models of thermospheric and ionospheric parameters to reproduce the observed morphology of the quiet-time zonal drifts for different seasons and distinct solar flux conditions. We focused, in particular, on an analytic representation (Haerendel et al. 1992) of the low-latitude electrodynamics that is commonly used to explain the observed behavior of the zonal plasma drifts (Eccles 1998;Chau and Woodman 2004;Rodrigues et al. 2012;Richmond et al. 2015).
We also focused on widely used empirical models of the thermospheric and ionospheric parameters to drive the analytic model. Noteworthy mentioning is that these models are part of the Pyglow python package that is readily available and commonly used by the geospace community. Pyglow includes the IGRF-12, IRI-2016, NRLMSISE-00, and three versions of the HWM (93, 07, and 14).
Finally, advances made in the ISR drift measurements at the Jicamarca Radio Observatory and long-term semiroutine operations (1986-2017) allowed observations of both vertical and zonal plasma drifts as a function of time and main F-region heights (240-560 km). While observations of the vertical drift were used to drive the analytic model, zonal drift observations were used to evaluate its performance for different seasons and distinct solar flux conditions. We were also able to evaluate the performance of the analytic representation for different wind model inputs.
Despite its simplicity, the analytic model can reproduce fairly well most of the features in the observed zonal plasma drifts, including the vertical shear of the zonal drifts associated with the evening plasma vortex.
During daytime hours the analytic model predicts similar results for the zonal drifts independently of the HWM used to drive the model. More importantly, the modeled drifts match well the behavior and magnitude of the observed daytime drifts for all seasons and solar flux conditions considered. This result indicates a good representation of the daytime neutral winds by the HWM models as well as a good representation of the low-latitude E-and F-region conductivities during the day.
During nighttime hours the analytic model results for different HWM inputs diverge, especially HWM07 and HWM14 from HWM93. Most of the results, however, are within the variability of the observations used to evaluate the models.
The nighttime results drive the overall performance of the analytic model, and we found that a single HWM cannot provide the best results for all seasons and solar flux conditions. The performance of the model outputs are quantified using a RMSE metric. HWM07 and HWM14 inputs tend to provide better results for LSF conditions, while HWM93 tend to perform better during HSF. These results can be explained, in most part, by a predominance of moderate-to-low solar flux observations used to create HWM07 and HWM14, and a solar flux parameterization implemented only in HWM93. The results indicate the necessity of additional work on modeling and validating the behavior of the nighttime winds and ionospheric conductivities. Additionally, following previous studies we neglected the contribution of the vertical current term (fourth term). This assumption has yet to be properly tested and requires the use of more complete ionospheric electric field model.
We also examined the main sources of zonal drift variability according to the analytic model (Eq. 6). Most of the morphology is controlled by the zonal wind dynamo term (first term). The contribution from meridional winds (second term) is negligible. Finally, the contribution from the vertical drifts (third term) cannot be neglected. This term makes a contribution to the daytime westward drifts on the order of the vertical drifts. It also contributes to the morphology of the zonal drifts in the bottomside F-region in the evening. The enhanced conductivity ratio at those heights acts as to weaken the shear in the zonal drifts.
We also examined the contribution from the E-and F-region to the zonal wind dynamo as predicted by the analytic model. The morphology of the zonal drifts in the region of observation (240-560 km altitude) is controlled in most part by the F-region winds. However, our results show that the E-region contribution is non-negligible during daytime hours, specially during December solstice. Also during December solstice, the E-region contribution to the F-region drifts is significant near sunset specially during LSF.
Finally, we would like to mention that the current work focused on evaluating the ability of current, readily available models of the thermosphere, ionosphere and electrodynamics to reproduce the observed zonal drifts. Given current limitations in the number of observations, we focused on evaluating the performance of the model for only two distinct solar flux conditions and specific DOYs representing each season. The encouraging results led to us to envision expanding the simple representation of the zonal drifts described here to cover other longitude sectors, and to accept DOY and a more complete range of solar flux conditions as inputs. The best approach and supporting data sets for such an effort are currently being investigated.

Appendix A: Collision frequency and conductivity
Collision frequencies were taken from Schunk and Nagy (2009). A recent study by Ieda (2020) provided updated coefficients for O + , NO + , and O + 2 collision frequencies. We have included the temperature-dependent resonant collisions in our calculations. The temperature independent ion-neutral collision frequencies are given in the form ν in = C in n n , where C in is the collision frequency coefficient for nonreasonant interactions, and the n n are the densities in cm −3 . The collision frequency coefficients are given in Table 2. Coefficients marked R represent resonant collisions. The temperature dependent resonant collisions are given by: where T r = (T i + T n )/2. Note that the definitions for the ν O+,O and ν O 2 +,O 2 resonant collision frequencies were taken from Ieda (2020). The final ion collision frequencies are then given by: The electron-neutral collision frequency is given by: where the summation is over all present ion species, e is the electric charge, B is the magnitude of the magnetic field, and n e and n i are the electron and ion densities in m −3 . The kappa term for species α is the ration of the gyro-to-collision frequency (κ α = ω α ν α = eB m α ν α ).