Equatorial spread F studies using SAMI 3 with two-dimensional and three-dimensional electrostatics

This letter presents a study of equatorial F region irregularities using the NRL SAMI3/ESF model, comparing results using a two-dimensional (2-D) and a threedimensional (3-D) electrostatic potential solution. For the 3-D potential solution, two cases are considered for parallel plasma transport: (1) transport based on the parallel ambipolar field, and (2) transport based on the parallel electric field. The results show that the growth rate of the generalized Rayleigh–Taylor instability is not affected by the choice of the potential solution. However, differences are observed in the structures of the irregularities between the 2-D and 3-D solutions. Additionally, the plasma velocity along the geomagnetic field computed using the full 3-D solution shows complex structures that are not captured by the simplified model. This points out that only the full 3-D model is able to fully capture the complex physics of the equatorial F region.


Introduction
Equatorial spread F (ESF) refers collectively to a family of plasma density irregularities that form in the equatorial F region ionosphere after sunset.The phenomenon has been extensively studied with the support of coherent scatter radars, ionosondes, airglow imagers, ground-based scintillation receivers, and instruments onboard rockets and satellites (see details in Woodman, 2009).
Topside ESF irregularities are driven by the generalized Rayleigh-Taylor instability near and above the F peak, and the plasma instability mechanism seems to be well under-stood.However, the first stages of ESF irregularity development, namely bottom-type and bottomside irregularities, and the plasma instabilities behind those two processes, have not been fully explained to date.Some of the mechanisms that have been proposed to explain the triggering of plasma waves are vertical shear in the horizontal plasma flow (Hysell and Kudeki, 2004), gravity waves (Huang and Kelley, 1996), and large-scale wave structure in the bottomside of the ionosphere (Tsunoda, 2005).
Numerical simulations have been used as a tool for the understanding of ESF.Three-dimensional models incorporating equipotential field lines (two-dimensional ionospheric potential solvers) have been extensively used (as, e.g., Huba et al., 2008 andRetterer, 2010).Two-dimensional ionospheric potential solvers use the fact that the plasma conductivity component parallel to the magnetic field is several orders of magnitude higher than the conductivity in the perpendicular direction, such that the magnetic field lines are assumed to be equipotential.This removes the dependence of the ionospheric potential in the dimension along the magnetic field, reducing the problem from 3-D to 2-D.Recent studies have shown that only 3-D models are able to capture the full nature of ESF (see, e.g., Aveiro and Hysell, 2010).Using both a 2-D and a 3-D ionospheric model, Aveiro and Hysell (2012) showed that the former is not able to simulate the initial stages of ESF and that only the 3-D model could simulate the 3 stages of ESF: bottom-type, bottomside, and topside.Their model was simplified in the sense that it only accounted for three ions: O + , O 2 + , and NO + .In this letter we describe two new upgrades to the Naval Research Laboratory (NRL) ionosphere model SAMI3/ESF: the use of a 3-D ionospheric electrostatic potential solver and the inclusion of a second-order transport scheme with flux limiters.We simulate the evolution of a topside ESF plasma irregularity using SAMI3/ESF with three different methods for the solution of the electrostatic potential: (1) 2-D potential (equipotential flux tubes), (2) 3-D potential with parallel transport using the ambipolar electric field associated with the parallel electron pressure (3-D-A), and (3) 3-D potential with parallel transport using the parallel electric field (3-D-F).The main goal is to evaluate how well the three approaches perform in terms of the physics captured and the complexity of the model for the simulation of plasma irregularities in ESF.

SAMI3/ESF model
The model was constructed using magnetic dipole coordinates (p, q, φ), where the tilt is matched to the magnetic declination in the longitude of interest.In our terminology, p represents the McIlwain parameter (L), q is the magnetic co-latitude, and φ is the longitude.The background current density including zonal wind forcing, plasma pressure, and gravity-driven currents is given as: where the terms on the right-hand side (RHS) represent the ohmic, diffusion, and gravity-driven currents, respectively.The terms ˆ , D, and ˆ represent the conductivity, diffusion, and gravity tensors, respectively (see, e.g., Shume et al., 2005 for an explicit definition of those terms).The variables U and B are the neutral wind speed and magnetic field, respectively.Perpendicular diffusion currents are neglected in the current version of the model, since their contribution is very small compared to the other two terms on the RHS.
The computation of the electrostatic potential is based on the solenoidal current density (∇ • J = 0) and can be written as where J = J 0 − ˆ • ∇ .In the 3-D model, Eq. ( 2) is solved directly in all 3 dimensions (p, q, φ).
In the 2-D model, the approach uses the fact that parallel conductivities are much larger than any of the perpendicular conductivities (Hall or Pedersen).Equation ( 2) is integrated along the parallel direction and the following condition is obtained: where the RHS is null, since parallel currents vanish at the bottom boundary of the ionosphere.This approach also indicates that there are no variations in the electrostatic potential along B, i.e., the field lines are equipotential.
The partial differential equation that describes the instantaneous behavior of the electrostatic potential based on the solenoidal current density for the 2-D approach then  simplifies to where and The terms p and φ are unit vectors, and σ P and σ H are the Pedersen and Hall conductivities, respectively.The h i coefficients are the scale factors that arise from the spherical to magnetic dipole coordinate system transformation (see, e.g., Huba et al., 2000), where the index i refers to the direction and represents the electrostatic potential.
Both Eqs. ( 4) and (2) are solved using the Bi-Conjugate Gradient Stabilized (BiCGStab) method using the SPARSKIT sparse matrix solver numerical library (Saad, 1990).The computational complexity increases with the matrix size and is O(N) for sparse matrix solvers.For instance, for a 3-D spatial grid n φ =96, n p =130, and n q =130 points, the 2-D potential matrix would be 12 480 × 12 480 and 1 622 400 × 1 622 400 for the full 3-D potential case.
Alongside the electrostatic potential, equations for the ion momentum and continuity are solved (see Huba et al., 2000 for a full description).The ion momentum equation in the parallel direction can be written as where V i , P i , n i , and m i represent the parallel velocity, plasma pressure, density, and mass for the i ion species, respectively.The terms ν in and ν ij are the collisions between the i ion species with the neutrals and the j ion species, respectively.The electric field component in the parallel direction comes from the electrostatic potential (E = −∇ ) in the 3-D-F model, and it is ambipolar (E = −∇ P e /n e ) in both the 2-D and 3-D-A models.Finally, temperature equations are not solved and instead thermodynamical equilibrium is assumed (T e = T i = T n ).A major simplification to the momentum equation in the perpendicular direction is made through the assumption that the main forcing is due to E ×B drifts.This condition will be relaxed in the next version of the model.The transport is performed using a Monotone Upwind Scheme for Conservation Laws (MUSCL) with the Sweby flux limiter (see, e.g., Trac and Pen, 2003) using a second-order Runge-Kutta method for the integration of the continuity equation in the perpendicular direction.The original SAMI3 model uses the donor cell method, which is only first-order accurate.

Model setup and results
The zonal boundaries of the model are periodic and the grid size is n φ = 96, n p = 130, and n q =130 points.The grid is ±2 • wide in longitude and spans between 90-690 km apex altitudes (L = 1.014-1.108).The simulation runs start at t 0 = 19:00 LT and the initial plasma conditions (densities and velocities) are obtained from a 36 h SAMI2 model run (Huba et al., 2000).A Gaussian perturbation is imposed on the plasma density at the bottomside of the F region.
Figure 1 depicts the initial conditions for the simulation runs at φ = 0.The left panels show diagnostics in the vertical direction at the magnetic equator.Zonal plasma drifts for all three models are similar below the F peak.The zonal shear node is located at ∼ 240 km, where the density gradient of the F region bottomside is very steep.At altitudes above ∼ 400 km the 2-D and the 3-D solutions differ, and the 2-D solution reflects the local zonal wind speeds (∼ 80 m s −1 ).Zonal plasma drifts for the 3-D models were slightly faster (∼ 90 m s −1 ), but were still slower than the maximum zonal wind speeds in the same flux tube (although not shown here, the magnitude of the zonal wind velocity reaches ∼ 120 m s −1 off-equator for these runs).This indicates that comparisons between local plasma drifts and neutral wind observations may be correlated, but differences of the order of 10 % might be expected (based on these simulation runs).Note that the higher in altitude the simulation goes, the larger the difference between the local neutral and plasma zonal velocities gets due to large off-equatorial Pedersen conductivities.
The right panels in Fig. 1 depict the parameters along the flux tube that intercept the Equator at 440 km altitude (L = 1.07).Since the variation of the electrostatic potential is the important parameter here, δ (= − • , where • is a constant arbitrary bias), is shown instead of .In the 3-D models, steep gradients in the potential along the flux tube were observed at ∼14 • magnetic latitude, where the low latitude E region was electrically decoupled from the equatorial F region below ∼ 150 km.Significant gradients in the electrostatic potential were observed as high as the F region bottomside.Since the parallel mobility is much larger (about 4 orders of magnitude) than the perpendicular mobilities at those altitudes, a small gradient in the potential along the magnetic field (i.e., a small parallel electric field) leads to a strong redistribution of plasma.The solution using the 2-D model assumes equipotential field lines such that the electrostatic potential does not vary along the flux tube.
Figure 2 shows the results at the end of the simulation (t ∼ 21:00 LT): (left) plasma densities in a cut through the equatorial plane and (right) O + ion parallel velocities in a cut through the meridional plane.From top to bottom, the solutions using the 2-D, 3-D with ambipolar electric field (3-D-A), and full 3-D model (3-D-F) runs, respectively.At 21:00 LT the instability was at the same stage of development in all model runs and the plasma irregularities had all crossed through the F peak (located at ∼ 380 km altitude).The plasma irregularity in the 2-D model was wider in longitude, indicating a zonal ballooning expansion.This indicates that the 2-D solution for the potential was smoother than in the 3-D models, and consequently the zonal plasma drifts were pointing outward from the center of the plasma irregularity over a larger zonal span.The opposite was observed in the 3-D models, where the more complex structure of the perpendicular electric fields and drifts was preserved, leading to short-scale forms and bifurcation (although not shown here, the plasma irregularities also bifurcated in the 2-D run, but at a later time step).The plasma density results for the 3-D models (3-D-A and 3-D-F) showed small differences above the F peak, but most of their structures remained similar.Efolding growth times estimated from the maximum vertical velocity time series were ∼ 14 min for all of the model runs, indicating that the generalized Rayleigh-Taylor instability is not affected by the choice of the potential solution.
Remarkable differences between the three models are seen in the O + ion parallel velocities in a cut through the meridional plane at the center of the plasma irregularity.Both 2-D and 3-D-A models presented a similar poleward plasma flow structure with steep gradients near the top of the plasma irregularity.The full model solution (3-D-F) showed additional structures in the parallel velocity that were not detected in the other two runs.Figure 3 shows a slice of the plasma densities and O + ion parallel velocities along the flux tube that intercepts the top of the plasma irregularity (at ∼440 km apex altitude).O + ion parallel velocities showed regions with enhanced equatorward flow at both the low latitude F valley and bottomside.Note that the simplified versions of the electrostatics model (2-D and 3-D-A) were not able to capture those features.The redistribution of the plasma density in the parallel direction for the 2-D model run was smoother than in the 3-D models, indicating enhanced parallel diffusivity.This seems to be caused by the absence of gradients in the perpendicular forcing along the magnetic field (∇ J 0⊥ ), since the 3-D-A model did not include any background current in the parallel direction (J 0 = 0), but was less diffusive than the 2-D model.

Discussions and conclusions
This letter presented a comparative study of the evolution of equatorial F region irregularities using a two-dimensional and a three-dimensional electrostatic potential solution in the NRL SAMI3/ESF model.For the 3-D potential solution, two cases are considered for parallel plasma transport: (1) transport based on the parallel electric field, and (2) transport based on the parallel ambipolar field.The results showed that the growth rate of plasma irregularities due to the Rayleigh-Taylor instability was unaffected by the choice of the ionospheric electrostatic potential approach.This was already expected, since the F region dynamo dominates during the nighttime, when the low latitude E region plasma conductivities are small.The growth rate of the generalized Rayleigh-Taylor instability is larger in collisionless regions and where the zonal background electric field is stronger, i.e., ∼ (g/ν in + E/B)n /n, and such parameters were not largely affected by the choice of the electrostatic potential approach.
However, the parallel flow was different.The typical poleward flow was observed when the potential solution did not include any forcing along the magnetic field, but a more structured flow with diverging and converging parallel flows was observed when the potential was fully solved in 3-D.Other differences included the lack of both fine structure and bifurcation in the 2-D solution when comparing the models at the same stage of evolution, which indicates that the 2-D solution for the ionospheric electrostatic potential does not fully describe the plasma flow structure of the ionosphere during equatorial spread F events, and therefore computing the full 3-D solution is necessary in order to capture the complex physics of the equatorial F region.

Fig. 2 .
Fig. 2. Plasma conditions at the end of the simulation (t ∼ 21:00 LT): (left) plasma densities in a cut through the equatorial plane and (right) O + ion parallel velocities in a cut through the meridional plane.From top to bottom, the solutions using the 2-D, 3-D-A, and 3-D-F models, respectively.

Fig. 3 .
Fig. 3. (Top) Plasma densities and (middle) O + ion parallel velocities along the flux tube that intercepts the top of the plasma irregularity (∼ 440 km apex altitude) for the (blue) 2-D, (green) 3-D-A, and (red) 3-D-F models.The bottom panel depicts the equivalent altitude (in km) at a given point.