Annales Geophysicae Three-dimensional simulation of equatorial spread-F with meridional wind effects

The NRL SAMI3 three-dimensional simulation code is used to examine the effect of meridional winds on the growth and suppression of equatorial spread F (ESF). The simulation geometry conforms to a dipole field geometry with field-line apex heights from 200 to 1600 km at the equator, but extends over only 4 degrees in longitude. The full SAMI3 ionosphere equations are included, providing ion dynamics both along and across the field. The potential is solved in two dimensions in the equatorial plane under a field-line equipotential approximation. By selectively including terms in the potential equation, the reduced growth predicted byMaruyama(1988) and the stabilization predicted byZalesak and Huba(1991) are separately realized. We find that ESF is stabilized by a sufficiently large constant meridional wind (60 m/s in our example).


Introduction
Equatorial spread F (ESF) was first identified by Booker and Wells (1938) as "diffuse echoes from the F-region of the ionosphere received continuously at night in equatorial regions over a wide range of wave-frequency" (Booker and Wells, 1938).Despite being first observed 70 years ago and the subject of intense research for the past 35 years (Haerendel, 1974;Ossakow, 1981;Kelley, 1989;Hysell, 2000), the underlying physical processes that control the onset and nonlinear evolution of ESF remain unknown.For example, the Correspondence to: J. Krall (jonathan.krall@nrl.navy.mil)cause(s) of day-to-day variability of ESF has not been established (Mendillo et al., 1992;Hysell and Kudeki, 2004;Tsunoda, 2005Tsunoda, , 2006)).In this light, ESF continues to be the most perplexing and fascinating ionospheric phenomenon at this time.Additionally, there is a societal impact attendant with ESF: the strong electron-density gradients associated with ESF affect the propagation of electromagnetic signals that can degrade global communication and navigation systems (de La Beaujardiere et al., 2004;Steenburgh et al., 2008).
The possibility that meridional winds could suppress ESF and contribute to the day-to-day variability of ESF has been raised by a number of researchers (Maruyama, 1988;Zalesak and Huba, 1991) based upon theoretical analyses and simplified numerical models.In particular, Maruyama (1988) computed the reduced growth rates that result both from redistribution of the Pedersen conductivity, along each field line, and from an overall increase in the Pedersen conductance associated with a meridional wind.In Maruyama (1988), the affect of the wind on the background plasma configuration in shown: a trans-equatorial meridional wind tends to push one peak in the Appleton anomaly upwards along the magnetic field and to push the conjugate peak downwards along the field.It is the affect of shifting this latter peak to lower altitude and higher neutral densities that increases both the local and the field-line-integrated Pedersen conductivity, reducing the computed ESF growth rate (Maruyama, 1988).In other words, Maruyama (1988) describes an indirect effect of the wind on the instability.Zalesak and Huba (1991) extended the analysis of Maruyama (1988) to consider the direct effect of the wind on the development of the instability (see also the Appendix below).In particular, a trans-equatorial wind directed along a fixed meridian provides wind components across the local magnetic field lines.In the instance where the Appleton anomaly peaks are shifted, as discussed above in reference to Maruyama (1988), the cross-field wind introduces an additional non-zero term into the potential equation.Zalesak and Huba (1991) showed that this additional effect can stabilize ESF.
Motivated by this work there have been a limited number of observational studies to determine the role of meridional winds on the occurrence of ESF, with mixed results.Mendillo et al. (1992) reported results of a two-night study of ESF using the ALTAIR radar and all-sky imaging at Kwajalein.They found that "surges" in the transequatorial winds on relatively short time scales (e.g., few hours to a day) could inhibit the onset of ESF.A subsequent study by Mendillo et al. (2001) during the Multi-Instrumented Studies of Equatorial Thermospheric Aeronomy (MISETA) campaign of 1998 found "no convincing evidence that meridional winds at dusk exert a strong influence on ESF onset."In a study of seasonal variability, Maruyama et al. (2009) find evidence for suppression of ESF by meridional winds.New wind and plasma data sets expected from the C/NOFS mission (de La Beaujardiere et al., 2004) are expected to shed light on this issue.
In this study, we apply a recently-developed threedimensional simulation model of ESF (Huba et al., 2008) to the issue of the partial or complete suppression of ESF by meridional winds.To our knowledge, this is the first comprehensive modeling study of the impact of meridional winds on the onset and evolution of ESF.

SAMI3 simulations
A modified version of the NRL 3-D global ionosphere code SAMI3 (Huba et al., 2008) is used in this analysis.SAMI3 models the plasma and chemical evolution of seven ion species (H + , He + , N + , O + , N + 2 , NO + and O + 2 ).The complete ion temperature equation is solved for three ion species (H + , He + and O + ) as well as the electron temperature equation.The plasma equations solved are as follows: The various coefficients and parameters are specified in Huba et al. (2000).Ion inertia is included in the ion momentum equation for motion along the geomagnetic field and E×B drifts are computed to obtain motion transverse to the field.
Neutral composition and temperature are specified using the empirical NRLMSISE00 model (Picone et al., 2002).The version of SAMI3 used here is modified relative to that used in past studies (Huba et al., 2005) in that the perpendicular electric field E ⊥ =−∇ is used self-consistently in SAMI3 to calculate the E×B drifts.SAMI3 uses a non-orthogonal, irregular grid that conforms to a model geomagnetic field.For this study the magnetic field is modeled as a dipole field aligned with the earth's spin axis, with the result that, in the equatorial plane, the SAMI3 grid is aligned with the coordinate system used in our analysis.For this study, the vertical and zonal neutral winds are set to zero and the meridional neutral wind is set to a constant value.Because the meridional wind crosses field lines, it provides nonzero values both parallel (V ns ) and perpendicular (V np ) to the field as shown in Fig. 1.
The potential equation is derived from current conservation (∇•J =0) in dipole coordinates (s, p, φ).The equation used in this study is where and Here g p is the component of g perpendicular to B, is the ion component of the Hall conductivity divided by the ion cyclotron frequency, is the Pedersen conductivity, pp = (p /b s )σ P ds, pφ = (1/pb s )σ P ds, =(1+3 cos 2 θ) 1/2 , i =eB/m i c, e =eB/m e c, B is the local geomagnetic field, B 0 is the geomagnetic field at the equator, b s =(r 3 E /r 3 ) , and r E is the radius of the earth.The field-line integrations above are along the entire flux-tube with the base of the field lines at 85 km.In deriving Eq. ( 5), we have neglected Hall-conductance terms relative to Pedersen-conductance terms and have neglected ∂/∂p terms relative to ∂/∂φ terms on the right-hand side.Equation ( 5) is similar to that derived in Haerendel et al. (1992).

Simulation parameters
The 3-D simulation model is initialized using results from the two-dimensional SAMI2 code.SAMI2 is run for 48 h using the following geophysical conditions: F10.7=170, F10.7A=170,A p =4, and day-of-year 263 (e.g., 20 September 2002).The geographic longitude is 0 • so universal time and local time are the same.The plasma is modeled from hemisphere to hemisphere up to ±26 • magnetic latitude; the peak altitude at the magnetic equator is ∼1600 km.The E×B drift in SAMI2 is prescribed by the Fejer/Scherliess model (Fejer and Scherliess, 1995).The plasma parameters (density, temperature, and velocity) at time 19:30 UT (of the second day) are used to initialize the 3-D model at each magnetic longitudinal plane.
The 3-D model uses a grid with magnetic apex heights from 200 km to 1600 km, and a longitudinal width of 4 • (e.g., 460 km) centered at 2 • .The grid is (n z , n f , n l )=(101, 130, 96) where n z is the number grid points along the magnetic field, n f is the number of "field lines," and n l is the number in longitude.This grid has a resolution of ∼10 km×5 km in altitude and longitude in the magnetic equatorial plane.The grid is periodic in longitude.In essence we are simulating a narrow "wedge" of the ionosphere in the post-sunset period.Finally, a Gaussian-like perturbation in the ion density is imposed at t=0: a peak ion density perturbation of 5% centered at 2 • longitude with a half-width of 0.25 • , and at an altitude z=400 km with a half-width of 60 km.The seed extends along the entire flux tube.

Meridional winds effects
As discussed in Maruyama (1988) and in the Appendix below, the indirect effect of the wind on ESF (Maruyama, 1988) enters through the first term on the right-hand side of Eq. ( 5) and the direct wind effect (Zalesak and Huba, 1991) enters through the second term, both of which are affected by integrations of wind-redistributed conductivities appearing on the left-hand side.Because a constant meridional wind crosses field lines (see Fig. 1), it provides nonzero values both parallel (V ns ) and perpendicular (V np ) to the field.Specifically, V ns acts through Eqs. ( 1) and (2), to redistribute the plasma, leading to interhemispheric asymmetries in the electron distribution and the Pedersen conductivity.
In this study, where the winds are northward, the southern ionosphere is pushed upwards and the conjugate northern ionosphere is pushed downwards.This is illustrated in Fig. 2 where contours of electron density are plotted in the latitudealtitude plane for zero wind (upper panel) and 60 m/s (lower panel).The corresponding effect on the Pedersen conductivity (see Fig. 3) is dramatic.As shown by Maruyama (1988), it is these asymmetries that affect the first term on the righthand side of Eq. ( 5) to reduce the growth rate.Below we refer to this as the "indirect" meridional wind effect.Zalesak and Huba (1991) showed that, in the presence of asymmetric n e and σ P distributions, such as shown in Figs. 2 and 3, V np acts directly through the ∂F φV /∂φ term in Eq. ( 5) to provide a net stabilizing contribution.For a constant meridional wind, this term is nonzero only if the interhemispheric plasma distribution is asymmetric.If the wind velocity is large enough, stability results (Zalesak and Huba, 1991).Because all ESF simulations necessarily include the ∂F φg /∂φ driving term, ESF simulations that also include the ∂F φV /∂φ term represent both the direct and indirect meridional wind effects.Below we refer to this as the "full" meridional wind effect.
As shown by Zalesak and Huba (1991) and in the Appendix, these meridional wind effects can be expressed in terms of field-line integrated growth rates: where and with L −1 n =∂ ln n 0 /∂p and Gravity being directed downwards, g p <0 and γ g is always positive (destablizing) in the bottomside F -layer.V np provides both positive and negative contributions (see Fig. 1) with the result that γ V is nonzero and negative (stabilizing) so long as the wind acts to shift n e in the direction of the wind as in Fig. 2.

Results
The initial condition for each simulation in this study is one in which the F-layer has been lifted to its highest point and is just beginning to fall.We performed several such runs of the SAMI3 code and present cases with constant meridional wind velocities of 0, 20, 40, 50, 60 and 80 m/s.In order to examine the indirect and full wind effects, we selectively exclude (indirect effect) or include (full wind effect) the second term on the right-hand side of Eq. ( 5).
We first consider the case with zero meridional winds.Figure 4 shows the growth of ESF over a 1.5 h period for no meridional winds.Both terms on the right-hand side of Eq. ( 5) are included, but have no stabilizing effect on the results.In order to quantify the growth of the instability, we tracked the bubble position, as indicated by the dots in Fig. 4. The dots indicate a position where the density is 80% of its local peak value.We also tracked the maximum upward (u p,max ) and horizontal (u φ,max ) components of the drift velocity, u≡E×B.
Figure 5 shows the evolution of the bubble height z (upper panel, long dashes).Figure 5 (middle panel) shows velocities dz/dt (dashed), u p,max (solid), and u φ,max (long dashes).As in Huba et al. (2008), we find peak upward velocities that are much larger than corresponding bubble velocities, indicating the generation of a substantial E field within the bubble.
Growth times (defined to be the e-folding time 1/γ ) can be estimated from the bubble velocity dz/dt or from either component of u.While each approach produces similar results, the vertical drift velocity u p,max gives the most consistent result during the early stage of the instability, where the growth is linear.Plots of log(dz/dt), log(u p,max ), and log(u φ,max ), are shown in Fig. 5 (lower panel).The growth time in this case, obtained via a linear fit to the first 40 min of the log(u p,max ) curve, is 24.1 min.This growth time is consistent with the 21.6 min growth time reported by Sultan (1996, see Fig. 3a) for no winds and a similar EUV parameter (sunspot number 140).As the instability develops nonlinearly, the growth rate increases and computed growth times decrease to match typical observed values of order 15 min.

Indirect meridional wind effect
Figure 6 shows the growth of ESF in a case with only the F φg term of Eq. ( 5) included and with a steady northward meridional wind of 60 m/s.As in Fig. 4, we tracked the bubble position, indicated by dots.Comparison of Figs. 6 and 4 show that the instability has slowed considerably, even during the nonlinear phase that is shown (in both cases the linear phase begins with a 5% perturbation at 19:30 UT).
Figure 7 shows the evolution of the bubble height z (upper panel) as well as the bubble velocity dz/dt and peak drift velocities u p,max and u φ,max (middle panel).As with Fig. 5, growth times can be estimated from plots of log(dz/dt), log(u p,max ), or log(u φ,max ), shown in the lower panel.In this case, the "bubble velocity" is initially subject to effects other than ESF, such as the evolution of the background F-layer.Based on a linear fit to the first 40 min of the log(u p,max ) curve, the growth time is 49.1 min.
As discussed above, the reduced ESF growth shown in Fig. 6 relative to that of Fig. 4 is caused by a reduction in the ESF driving term ∂F φg /∂φ.This is shown in Fig. 8, which shows plots of ∂F φg /∂φ at altitude 415 km, where the driving term was strongest during the linear growth phase of the instability.The driving term is plotted at times 19:30 UT (dashed), 20:15 UT (long dashes) and 21:00 UT (solid) for the zero wind (upper panel) and 60 m/s (lower panel; scale reduced by a factor of 10) "indirect wind effect" cases.While the driving term is clearly much smaller in the 60 m/s case relative to the zero wind case, it does grow with time.The F-layer remains unstable, albeit with a longer growth time.
Results for these "indirect wind effect" cases (with only the F φv term included in Eq. 5) are summarized in Fig. 9, which shows line plots of log 10 (u p,max ) versus local time for no wind (solid), 20 m/s (dashes), 40 m/s (long dashes), 60 m/s (dash-dot), and 80 m/s (lower solid line), with velocities in units of m/s.Clearly the growth of the instability is reduced with rising wind speed, consistent with past results (Maruyama, 1988;Sultan, 1996).

Full meridional wind effect
Similar to the indirect-effect runs discussed above, we performed a number of SAMI3 ESF full-wind-effect simulations with both the F φg and F φV terms included in Eq. ( 5) Fig. 10.Line plots showing log 10 (u p,max ) versus local time for no wind (solid), 20 m/s (dashes), 40 m/s (long dashes), and 50 m/s (dash-dot), and 60 m/s (lower solid line) for the "full wind effect" runs.Wind velocities differ from those used in Fig. 9. and with constant meridional wind velocities of 0, 20, 40, 50 and 60 m/s.Results are summarized in Fig. 10, which shows log 10 (v max ) versus time for each case.Comparing Figs. 9 and 10 shows that, as predicted (Zalesak and Huba, 1991;Sultan, 1996), the instability is stabilized with a sufficiently large meridional wind.In the 60 m/s case, velocities associated with the initial perturbation simply die down and a quiescent state is obtained.
The driving terms ∂F φg /∂φ (dashed) and ∂F φV /∂φ (long dashes) are plotted in Fig. 11 for the 60 m/s case, in which ESF was completely stabilized.These are plotted at an altitude of 415 km and at times 19:30 UT (lower panel), 21:00 UT (middle panel), and 22:30 UT (upper panel).This figure shows that, while the individual driving terms increase with time, ∂F φV /∂φ opposes ∂F φg /∂φ in such a way that the net driving term (solid curve) decreases versus time.

Summary
Figure 12 shows that ESF growth times are increased via the indirect wind effect (solid squares) and further increased with the inclusion of the direct wind effect (dots).These results show that the reduced growth predicted by Maruyama (1988) and the stabilization predicted by Zalesak and Huba (1991) are separately realized and that a constant meridional wind of 60 m/s stabilizes spread-F.For comparison to these results, we have computed growth times based on the approximate analytic formulas, Eqs.(10-12).These growth times, shown as open symbols in Fig. 12, are computed for each simulation during ESF linear growth (one hour after each simulation be-Fig.11.ESF driving terms ∂F φg /∂φ (dashed) and ∂F φV /∂φ (long dashes) at altitude 415 km and at times 19:30 UT, 21:00 UT, and 22:30 UT for the 60 m/s "full wind effect" case.The net driving term is the solid curve in each panel.gins).We see that the analytic formulas produce reasonable results for low and moderate wind speeds.Above 40 m/s, analytic growth times do not vary as strongly with wind speed as was seen in the simulations.For completeness we note that growth times computed using Eqs.(10-12) vary as ESF develops.Typically, a computed growth time exhibits variations of order 10% while the corresponding simulation exhibits constant exponential growth.
These results are consistent with those of Sultan (1996).To be clear, what is reported by Sultan (1996) are computed growth rates based on a parameterized ionospheric density model, FAIM (fully analytic ionospheric model), that has empirically-determined winds as one of its inputs instead of the constant wind used in our study.Thus, an appropriate comparison is between cases in which the SAMI3 and FAIM wind-affected n e distributions are similar as is the case with Fig. 5b of Sultan (1996), where the growth time is 59.5 min with all wind effects included, and our 40 m/s full-windeffect case (latitude-altitude distribution not shown) in which the growth time is 52.3 min.Even though the n e distribution shown in Fig. 5b of Sultan (1996) is similar to our 40 m/s case, the Sultan (1996) result features a much larger peak local wind velocity of 200 m/s.However, it is the asymmetry in the plasma distribution, and the corresponding asymmetry in the Pedersen conductivity, that produces the indirect wind effect (reducing growth) and enables the direct wind effect (further reducing growth and leading to possible stabilization).A very large meridional wind is of no consequence without these inter-hemispherical asymmetries.The 60 m/s wind case shown in Figs. 2 and 3 tilts the electron distribution northward to a degree that is clearly greater than that shown in Sultan (1996, Fig. 5b).In other words, our steady 60 m/s wind shows a stronger effect on both the longitudealtitude n e distribution and on ESF growth than was found for more realistic wind pattern in which the peak velocity is 200 m/s.
In order to focus on leading-order meridional-wind effects, a number of other phenomena have neglected in this study, such as those associated with zonal winds and those associated with second-order terms in Eq. ( 5).Specifically Hall-conductance terms were neglected relative to Pedersenconductance terms and ∂/∂p terms were neglected relative to ∂/∂φ terms on the right-hand side of Eq. ( 5).Additional simulations show that the baseline zero-wind case is only minimally affected when the Hall-conductance and ∂/∂p terms are included in the potential equation.
A number of interesting questions have been left for future study, such as the effect of zonal winds on ESF dynamics, the effect of more realistic wind patterns, the determination of the maximum bubble-height for various conditions, and the three-dimensional shape and dynamics of the result- ing high-altitude depletions.We also have not considered more complex initial perturbations, such as those that lead to bifurcations (Huba et al., 2008), or a higher-order transport scheme that would better capture such physics (see e.g.Huba and Joyce, 2007).
In particular, we have been careful to use similar initial conditions in all cases, with the same specified E×B drift (Fejer and Scherliess, 1995) being used in each case to lift the F-layer prior to the initiation of the self-consistent SAMI3 simulation and with each simulation being initiated at nearly the same local time.Thus, the importance of the E×B drift magnitude on these results has not been systematically evaluated.A study by Mendillo et al. (2001) strongly suggests that ESF growth varies more strongly with the strength of the eastward electric field than with the strength of the meridional wind.
We did find that ESF growth is sensitive to the initial state of the F-layer as determined by the ESF initiation time in each simulation.By this, we mean the time at which the imposed drift function used in the SAMI2 initialization is replaced by the self-consistently computed potential and resulting drifts of the SAMI3 simulation.For example, ESF simulations initiated at 19:00 local time (equatorial F-layer lifted up to 350 km before initiation of ESF in the zero-wind case) grew about half as quickly as those initiated at 19:30 (F-layer lifted up to 400 km).

Conclusions
We have used the SAMI3 three-dimensional simulation code to examine the effect of meridional winds on the growth and suppression of equatorial spread-F.The simulation geometry conforms to a dipole field geometry that extends from 200-1600 km at the equator, down to and altitude of 85 km at the ends of each field line, and over 4 • in longitude (longitudinal width 460 km).The full SAMI3 ionosphere equations are included, providing ion dynamics both along and across the

Fig. 1 .
Fig. 1.Diagram indicating field lines, a uniform northward meridional wind, and wind components along and across the field.

Fig. 2 .
Fig. 2. Contours of electron density plotted in the latitude-altitude plane for zero wind (upper panel) and 60 m/s (lower panel).These show initial states, before ESF develops.

Fig. 3 .
Fig. 3. Contours of Pedersen conductivity plotted in the latitudealtitude plane for zero wind (upper panel) and 60 m/s (lower panel).These show initial states, before ESF develops.

Fig. 4 .
Fig. 4. Contours of n e showing ESF development with no meridional wind.

Fig. 6 .
Fig. 6.Contours of n e for the 60 m/s "indirect wind effect" case.