SAMI3_ICON: Model of the Ionosphere/Plasmasphere System

Abstract

The NRL ionosphere/plasmasphere model SAMI3 has been modified to support the NASA ICON mission. Specifically, SAMI3_ICON has been modified to import the thermospheric composition, temperature, and winds from TIEGCM-ICON and the high-latitude potential from AMIE data. The codes will be run on a daily basis during the ICON mission to provide ionosphere and thermosphere properties to the science community. SAMI3_ICON will provide ionospheric and plasmaspheric parameters such as the electron and ion densities, temperatures, and velocities, as well as the total electron content (TEC), peak ionospheric electron density (NmF2) and height of the $F$ layer at NmF2 (hmF2).

1 Introduction

The ICON mission is designed to provide coverage of the ionosphere/thermosphere system in the low- to mid-latitude region. The ICON satellite will have a circular orbit at 575 km with a 27^∘ inclination. One of the primary goals of the mission is to measure the thermospheric winds in the $E$ and $F$ regions of the ionosphere. The nominal latitudinal range of wind measurements is 12^∘ S to 45^∘ N. These winds are critical to modeling the dynamics of the ionosphere because of ion-neutral coupling and the generation of the low- to mid-latitude electric field driven by the neutral wind dynamo. In support of this mission, the NRL ionosphere/plasmasphere model SAMI3 has been modified to incorporate thermospheric composition, temperature, and winds from the TIEGCM-ICON model that will ingest data from the ICON satellite. More information about TIEGCM can be found in Qian et al. (2014) and Richmond and Maute (2014) (and references therein). Additionally, SAMI3 will use the high-latitude potential from the AMIE model (Richmond et al. 1992; Crowley and Hackert 2001). In this paper we present a brief description of the SAMI3_ICON model and preliminary results for the May 9, 2005 magnetic storm.

2 Model Description

SAMI3 (Sami3 is A Model of the Ionosphere) is a three-dimensional, physics-based model of the ionosphere developed at the Naval Research Laboratory (NRL). It is based on SAMI2 (Huba et al. 2000) a two-dimensional model of the ionosphere. SAMI3 models the plasma and chemical evolution of seven ion species (H⁺, He⁺, N⁺, O⁺, $\mathrm{N}_{2}^{+}$, NO⁺ and $\mathrm{O}_{2}^{+}$). The temperature equation is solved for three ion species (H⁺, He⁺ and O⁺) and for the electrons. Ion inertia is included in the ion momentum equation for motion along the geomagnetic field and ${\mathbf{E}} \times {\mathbf{B}}$ drifts are computed to obtain motion transverse to the field. SAMI3 includes 21 chemical reactions and radiative recombination.

The basic plasma equations solved by SAMI3 are the following.

• 1. Ion continuity:

  $$ \frac{\partial n_{i}}{\partial t}+\nabla \cdot ( n_{i}\mathbf{V} _{i} ) = \mathcal{P}_{i}-\mathcal{L}_{i}n_{i} $$

  (1)

• 2. Ion velocity:

  $$\begin{aligned} \frac{\partial \mathbf{V}_{i}}{\partial t}+ \mathbf{V}_{i}\cdot \nabla \mathbf{V}_{i}&= -\frac{1}{\rho_{i}}\nabla \mathbf{P}_{i}+ \frac{e}{m _{i}}\mathbf{E}+ \frac{e}{m_{i}c} \mathbf{V}_{i}\times \mathbf{B+g} \\ &\quad{}-{\nu_{in}}(\mathbf{V}_{i}-\mathbf{{V}_{n}}) -\sum_{j}\nu_{ij} ( \mathbf{V}_{i}- \mathbf{V}_{j} ) -\nu_{a}(V_{i\parallel }) V_{i \parallel } \end{aligned}$$

  (2)

• 3. Ion temperature:

  $$ \frac{\partial T_{i}}{\partial t}+\mathbf{V}_{i}\cdot \nabla T_{i}+ \frac{2}{3}T_{i}\nabla \cdot \mathbf{V}_{i} + \frac{2}{3}\frac{1}{n _{i}k}\nabla \cdot \mathbf{Q}_{i}=Q_{in}+Q_{ij}+Q_{ie} $$

  (3)

• 4. Electron velocity along $\mathbf{B}$:

  $$ 0=-\frac{1}{n_{e}m_{e}}b_{s}\frac{\partial P_{e}}{\partial s}-\frac{e}{m _{e}}E_{s} $$

  (4)

• 5. Electron temperature:

  $$ \frac{\partial T_{e}}{\partial t}-\frac{2}{3}\frac{1}{n_{e}k}b_{s} ^{2}\frac{\partial }{\partial s}\kappa_{e}\frac{\partial T_{e}}{ \partial s}=Q_{en}+Q_{ei} + Q_{phe} $$

  (5)

In the above, $n_{i}$ is the ion density, $\mathbf{V}_{i}$ is the ion velocity, $\mathbf{V}_{n}$ is the neutral wind velocity, $T_{i,e}$ is the ion/electron temperature, $\mathcal{P}_{i}$ is the ion production rate (i.e., photoionization and chemistry), $\mathcal{L}_{i}$ is the ion loss rate (i.e., chemistry). $P_{i.e}$ is the ion/electron pressure, $\mathbf{E}$ is the electric field, $\mathbf{B}$ is the geomagnetic field, $\mathbf{g}$ is gravity, $\nu_{in}$ is the ion-neutral collision frequency, $\nu_{ij}$ is the ion-ion collision frequency, $\nu_{a}$ is a velocity-dependent collision term that prevents ions from becoming supersonic along the magnetic field, $k$ is the Boltzmann constant, $Q_{i}$ is the ion heat flux, $\kappa_{e}$ is the electron thermal conductivity, $s$ is the coordinate along $\mathbf{B}$, $b_{s} = B_{s}/B _{0}$ where $B_{0}$ is the earth’s surface magnetic field, $Q_{\alpha ,\beta }$ is the heat exchange rate between species, $\alpha$ and $\beta$, and $Q_{phe}$ is the photoelectron heating rate. The explicit definitions of parameters (e.g., $\nu_{in}$, $P_{i}$, $L_{i}$, etc.) are defined in Huba et al. (2000).

The neutral composition, temperature, and winds can be specified in SAMI3 by analytical models, empirical models (NRLMSISE00 (Picone et al. 2002), HWM93 (Hedin et al. 1991), HWM14 (Drob et al. 2015)) or physics-based models (e.g., TIEGCM, TIMEGCM). In SAMI3_ICON, we use thermospheric inputs from TIEGCM-ICON. At altitudes below the lower TIEGCM-ICON boundary (i.e., 97 km) we use NRLMSISE00 and HWM14 to specify the thermosphere. We have tested the code by varying the lower boundary of SAMI3 below and above 97 km, and found that the transition to NRLMSISE00 and HWM14 does not affect the ionospheric dynamics. At altitudes above the upper TIEGCM-ICON boundary (600–700 km) we logarithmically extrapolate the neutral density, and use a constant temperature and neutral wind. Additionally we use NRLMSISE00 to specify the neutral hydrogen density which is not provided by TIEGCM-ICON.

SAMI3_ICON uses a nonorthogonal, nonuniform, fixed grid with one axis aligned with the geomagnetic field that is ‘seamless’ (i.e., there are no transition regions between low- and high-latitude). In SAMI3_ICON we use the Richmond apex model (Richmond 1995) consistent with that TIEGCM-ICON and AMIE. The latitudinal extent (magnetic) of SAMI3_ICON can be up to $\pm 89^{\circ }$. The grid is designed to optimize the numerical mesh so that the spatial resolution decreases with increasing altitude. The plasma is transported transverse to the geomagnetic field using a finite volume method in conjunction with the donor cell method. The code is parallelized in the zonal direction using the Message Passing Interface (MPI) method which is critical to computational efficiency for large-scale simulations.

The electric field in SAMI3_ICON is found by solving the potential equation in dipole coordinates ($q,p,\phi$) in the magnetic equatorial plane ($p,\phi$):

$$\begin{aligned} &{\underbrace{ { \frac{\partial }{\partial p}p\varSigma_{pp} \frac{\partial \varPhi }{\partial p} +\frac{\partial }{\partial \phi }\frac{ \varSigma_{p\phi }}{p} \frac{\partial \varPhi }{\partial \phi } } }_{ {\mathrm{Pedersen}}} } \quad {\underbrace{ { -\frac{ \partial }{\partial p} \varSigma_{H} \frac{\partial \varPhi }{\partial \phi } +\frac{\partial }{\partial \phi }\varSigma_{H} \frac{\partial \varPhi }{\partial p} } }_{ {\mathrm{Hall}}} } \\ &\quad= {\underbrace{ { \frac{\partial F_{pV}}{\partial p} +\frac{ \partial F_{\phi V}}{\partial \phi } } }_{ {\mathrm{neutral \ winds}}} } \quad {\underbrace{- \frac{\partial F_{pg}}{\partial p} +\frac{ \partial F_{\phi g}}{\partial \phi }} _{ {\mathrm{gravity}}}} \end{aligned}$$

(6)

where the parameters are defined in the Appendix (Huba and Joyce 2010). This equation is based on two assumptions: $\nabla \cdot {\mathbf{J}} = 0$ and equipotential field lines. In SAMI3_ICON, the above equation is solved using the neutral wind specified by TIEGCM-ICON, as well as the neutral composition and temperature, to calculate the dynamo electric field. The Richmond apex model is used to defined the thermospheric parameters at the appropriate geographic location and geometric factors relative to the dipole magnetic field used in the solution of Eq. (6).

The high latitude electric field is determined from the AMIE potential. The potential used in the model is the sum of these two potentials. There is a minor issue in that the AMIE potential implicitly contains a contribution from the neutral wind dynamo that SAMI3 calculates. However, the neutral wind dynamo contribution to the potential at high latitude is generally very small and does not significantly affect ionospheric dynamics.

It is worth noting that SAMI3 can specify a non-uniform longitudinal grid so that over a localized longitudinal region (e.g., 60^∘) we can achieve a grid resolution of 0.0625^∘ ($\sim 7$ km) to self-consistently study the development of equatorial plasma bubbles (Huba and Joyce 2010). We know of no other global ionosphere model with this capability.

TIEGCM is a 3D coupled thermosphere-ionosphere electrodynamic general circulation model. It is based on TI-GCM that was developed by Roble et al. (1988). TI-GCM was extended to include self-consistent electrodynamic interactions between the ionosphere and thermosphere (Richmond et al. 1992). Although not used here, it is worth mentioning that another version of the model has been extended down to 30 km, to include the mesosphere and upper stratosphere and is known as TIMEGCM (Roble and Ridley 1994). TIEGCM predicts the thermospheric winds, temperatures, major and minor composition, electron densities and electrodynamic quantities globally from 95 km to about 500 km altitude (altitude of top pressure level depends on solar cycle). The current TIEGCM model uses a fixed geographic grid with a 2.5^∘ latitude $\times\ 2.5^{\circ }$ longitude horizontal resolution, and a vertical resolution of a quarter pressure scale height. The model time-step is typically 3 minutes, but is reduced to 30 sec or less in order to model more active conditions. The inputs required by the TIEGCM include the solar flux (parameterized by the F10.7 flux), auroral particle precipitation, high latitude electric fields, and tides propagating up from below the 95 km lower boundary. Various authors have been successful using the TIEGCM in simulating thermospheric composition and density changes caused by spatially and temporally varying energy and momentum inputs into the high latitude thermosphere-ionosphere system (Crowley et al. 2008; Crowley et al. 2010). The new feature of TIEGCM-ICON to be used for the ICON mission is that the lower boundary will be modulated by Hough mode extensions obtained from analysis of ICON wind data. This will provide more accurate inputs for the complex tidal motions affecting thermospheric winds.

In general, all of the desired electrodynamic quantities that are required to drive the TIEGCM-ICON or SAMI3_ICON model are not available globally and simultaneously. Therefore, we generally drive the TIEGCM-ICON with the time-dependent high-latitude ionospheric electric potential and auroral precipitation patterns provided by the Assimilative Mapping of Ionospheric Electrodynamics (AMIE) procedure (Richmond et al. 1992; Crowley and Hackert 2001). AMIE is an inversion technique that ingests data from a diverse range of sources to produce a realistic representation of the high-latitude electrodynamic state for a given time. The data inputs typically include electric fields derived from ion velocities measured by radars, satellites, and digital ionosondes, together with magnetic perturbations from ground- and space-based magnetometers including AMPERE (Wilder et al. 2012). Using these data, the distribution of various electrodynamic parameters such as the electric potential and electric field, height-integrated conductivity, and horizontal and FACs can be derived through the electrodynamic equations. Thus, combining AMIE with TIEGCM-ICON and SAMI3_ICON provides an extended capability to perform detailed investigations that obtain better agreement between the observed and modeled atmospheric responses (Crowley et al. 2010), thereby helping to advance the modeling component of the ICON mission.

3 Preliminary Results

We show results from SAMI3_ICON for the time period May 7-9, 2005. For these days SAMI3_ICON reads in TIEGCM-ICON netcdf data files for thermosphere composition, temperature, and winds as well as AMIE potential netcdf data files. We have chosen a stormtime period to capture dynamic effects associated with the AMIE potential absent during quiet times. The simulations are done for each day separately (i.e., for 24 hrs) to test the restart capability in SAMI3_ICON. Here, for example, we initialize the May 8 simulation at 00:00 UT with the simulation results from the May 7 simulation at 24:00 UT. This ensures continuity of simulation results from day-to-day.

The geophysical inputs to SAMI3_ICON are $\mbox{F10.7} = 111$, $\mbox{F10.7A} = 96$, and $\mbox{Ap} = 4$. The grid is $(nz,nf,nl) = (204,80,96)$ where $nz$ is the number of grid points along the magnetic field, $nf$ is the number of ‘field lines’ (i.e., the latitudinal resolution in each hemisphere; there are $nf$ magnetic field footpoints in each hemisphere), and $nl$ is the number of grid points in the longitude. The resolution along and across the magnetic field is variable, but is $\lesssim 20$ km in the $F$ region.

In Fig. 1, the solar wind magnetic field component $B_{z}$ and the Dst are shown as a function of time (UT) for the period of the simulation study. Although the IMF $B_{z}$ experienced large fluctuations, there were two periods of southward $B_{z}$ on average: ∼ 16:00 (day 127) to ∼ 12:00 (day 128) (storm 1) and ∼ 12:00 (day 128) to ∼ 20:00 (day 128) (storm 2). These two time periods led to modest magnetic storm activity as seen in the Dst with two main phases: ∼ 20:00 (day 127) to ∼ 02:00 (day 128) and ∼ 12:00 (day 128) to ∼ 20:00 (day 128).

Fig. 1

The solar wind magnetic field $B_{z}$ and Dst as a function of time (UT) for May 7–9, 2005 (days 127, 128, 129)

In Figs. 2–6 we show contour plots of log₁₀(TECU) for both GPS data (top panels) and SAMI3_ICON results (bottom panels). The contour plots use the same color scale and roughly the same longitude/latitude mapping in the low- to mid-latitude sector (i.e., within $\pm 60^{\circ }$). The white contour lines in the SAMI3_ICON plots are the electrostatic potential from the neutral wind dynamo and AMIE. AMIE potentials are computed independently for both the northern and southern hemispheres, but in this study the AMIE potential in the northern hemisphere is mapped to the southern hemisphere. SAMI3_ICON will be upgraded to allow for different AMIE potentials in northern and southern polar regions. The plots are spaced in 12 hr intervals starting on March 7, 2015 (day 127) at 18:00 UT (Fig. 2) and ending on March 9, 2015 (day 129) at 18:00 UT (Fig. 6). (Note: The white circle in the lower right portion of the SAMI3_ICON plots is the magnetic south pole; there is no data here because the model only extends to $\pm 88^{\circ }$ in this study.)

Fig. 2

Global plot of $\log_{10}(\mathrm{TECU})$ for May 7, 2005 (day 127) at 18:00 UT: GPS data (top) and SAMI3_ICON results (bottom). The potential contours are in the range −95 kV to 49 kV

Figure 2 is at 18:00 UT on March 7, 2015 (day 127) and is at the beginning of the sudden commencement of the magnetic storm. The region of largest TEC is around the magnetic equator in the South American and mid-Atlantic sectors. The model results agree quite well with the data in these regions. The GPS data has somewhat higher TEC in the North American sector than the model, while the model has somewhat higher TEC in the Indian and Australian sectors. These differences may be due to a difference between the electric field used in the model versus the actual electric field, as well as differences in the thermospheric composition and winds. The electrostatic potential is largest in the high-latitude region as expected; there is a weak coupling of the mid-latitude potential to the low-latitude as indicated by the contour line penetrating to the equator.

Figure 3 is at 06:00 UT on March 8, 2015 (day 128) and is in the beginning of the recovery phase of the ‘first’ magnetic storm. Overall there is good agreement between the data and model. Interestingly, in the data there is a latitudinally narrow, low TEC region extending across lower Canada that is also apparent in the simulation results.

Fig. 3

Global plot of $\log_{10}(\mathrm{TECU})$ for May 8, 2005 (day 128) at 06:00 UT: GPS data (top) and SAMI3_ICON results (bottom). The potential contours are in the range −75 kV to 68 kV

Figure 4 is at 18:00 UT on March 8, 2015 (day 128) and is at the end of the main phase of the ‘second’ magnetic storm. Again, there is good agreement between the data and model. In the data there is a of plume extending from the northeast United States into northern, central Canada that is also apparent in the simulation results. This plume is caused by the convection of plasma from mid-latitudes into the polar cap under the influence of the AMIE potential which extends to lower latitudes during the storm. This demonstrates the value of including high-fidelity potential from AMIE in the simulation code.

Fig. 4

Global plot of $\log_{10}(\mathrm{TECU})$ for May 8, 2005 (day 128) at 18:00 UT: GPS data (top) and SAMI3_ICON results (bottom). The potential contours are in the range −75 kV to 68 kV

Figure 5 is at 06:00 UT on March 9, 2015 (day 129) and is well into the recovery phase of the ‘second’ magnetic storm. There is good agreement between the data and model especially in the Indian sector. As in Fig. 3, there is a latitudinally narrow, low TEC region extending across lower Canada but is more intense than that in Fig. 3. This feature is also apparent in the simulation results and is indicative of the mid-latitude trough(Rodger 2008). Additionally, there is a region of low TEC extending from the west coast of Africa to the eastern coast of South America below the equator. This behavior is also seen in the simulation results.

Fig. 5

Global plot of $\log_{10}(\mathrm{TECU})$ for May 9, 2005 (day 129) at 06:00 UT: GPS data (top) and SAMI3_ICON results (bottom). The potential contours are in the range −59 kV to 32 kV

Figure 6 is at 18:00 UT on March 9, 2015 (day 129) and is further into the recovery phase of the ‘second’ magnetic storm. The ionosphere is evolving into a quiet time state and there is very good agreement between the data and model in the low- to mid-latitude ionosphere where ICON will be making measurements.

Fig. 6

Global plot of $\log_{10}(\mathrm{TECU})$ for May 9, 2005 (day 129) at 18:00 UT: GPS data (top) and SAMI3_ICON results (bottom). The potential contours are in the range −59 kV to 32 kV.

Finally, In Fig. 7 we show isocontours of the electron density in the plasmasphere at times (1) prior to the magnetic storm at 18:00 UT on day 127 and (2) in the recovery phase of the storm at 18:00 UT on day 129. The green contour is at $n_{e} = 10^{2}\mbox{ cm}^{-3}$ and the orange contour is at $n_{e} = 10^{6}\mbox{ cm}^{-3}$. In the left panels the view is from the sun, while in the right panels the sun is at the bottom. In the pre-storm period (top panels) the plasmasphere extends to ∼ 6 $R_{E}$ in the equatorial plane (top right panel) and has a toroidal shape. In the recovery phase period (bottom panels) the plasmasphere has undergone stormtime erosion and now extends to only $\sim 4 R_{E}$. Moreover, the plasmasphere becomes distorted in the equatorial plane and there is a bulge in the afternoon sector (bottom right panel). These results are consistent with observations in general, although we do not see a clear plasmaspheric plume develop in this simulation study (Goldstein et al. 2005; Huba and Sazykin 2014). Recently, Krall et al. (2014) presented results that suggest thermospheric winds have an impact on the shape and structure of the plasmasphere during quiet times. Thus, we should be able to assess this hypothesis with ICON wind data and plasmasphere observations.

Fig. 7

Isocontours of the electron density in the plasmasphere at times 18:00 UT (day 127) and 18:00 UT (day 129). The green contour is at $n_{e} = 10^{2}\mbox{ cm}^{-3}$ and the orange contour is at $n_{e} = 10^{6}\mbox{ cm}^{-3}$. In the left panels the view is from the sun, while in the right panels the sun is at the bottom

4 Summary

In this paper we have presented a brief description of the SAMI3_ICON model and preliminary results for the May 9, 2005 magnetic storm. The simulation results for this storm are in reasonably good agreement with GPS TEC data, and the overall behavior of the plasmasphere is consistent with observations. The key aspect of this model is that we use input data from TIEGCM-ICON and AMIE via netcdf files for the neutral composition, temperature, winds, and high-latitude potential. These models will be run on a daily basis during the ICON mission using, in part, data from the ICON satellite. One of the main reasons that the model results do not agree with observational data in the low- to mid-latitude ionosphere is the specification of the thermospheric neutral wind. A primary purpose of the ICON mission is to measure the the neutral wind on a global scale in the low- to mid-latitude $E$ and $F$ regions. Incorporating the measured winds into ionospheric models, such as SAMI3_ICON, is expected to improve the agreement between models and data.

SAMI3_ICON data will be post-processed on a daily basis and output to a common ICON format for netcdf files; the data that will be output will be the global 3D electron density, ion densities, velocities, and temperatures, as well as the TEC, NmF2, and hmF2. These data products (Level 4) will be available to the scientific community on the ICON mission website on a daily basis.

Appendix

The definitions of the terms in Eq. (6) are as follows in dipole coordinates where $R_{E}$ is the radius of the earth.

$$\begin{gathered} q = \frac{r_{E}^{2}{\mathrm{cos}}\theta }{r^{2}} \qquad p = \frac{r}{r_{E} {\mathrm{sin}}^{2}\theta } \qquad \phi = \phi \\ \varSigma_{pp} = \int (\sigma_{P}/b_{s})\Delta \, ds \qquad \varSigma_{p \phi }= \int (\sigma_{P}/b_{s}\Delta )\, ds \qquad \varSigma_{H} = \int ( \sigma_{H}/b_{s})\, ds \\ F_{pV} = \int \frac{B_{0}}{c} r\sin \theta (\sigma_{P}V_{n\phi }+ \sigma_{H}V_{np})ds \\ F_{pg} = \int r\sin \theta \frac{B_{0}}{c}\frac{1}{\varOmega_{i}}\sigma _{Pi}g_{p} ds \\ F_{\phi V} = \int \frac{B_{0}}{c} \frac{r_{E}\sin^{3}\theta }{\Delta } (-\sigma_{P}V_{np} + \sigma_{H}V_{n\phi }) ds \\ F_{\phi g} = \int \frac{r_{E}\sin^{3}\theta }{\Delta }\frac{B_{0}}{c} \frac{1}{ \varOmega_{i}} \sigma_{Hi}g_{p} ds \\ \sigma_{P} = \sum_{i} \frac{n_{i}ec}{B}\frac{\nu_{in}/\varOmega_{i}}{1 + \nu_{in}^{2}/\varOmega_{i}^{2}} + \frac{n_{e}ec}{B}\frac{\nu_{en}/\varOmega _{e}}{1 + \nu_{en}^{2}/\varOmega_{e}^{2}} \\ \sigma_{H} = -\sum_{i} \frac{n_{i}ec}{B}\frac{1}{1 + \nu_{in}^{2}/ \varOmega_{i}^{2}} + \frac{n_{e}ec}{B}\frac{1}{1 + \nu_{en}^{2}/\varOmega _{e}^{2}} \\ \sigma_{Hi} = -\sum_{i} \frac{n_{i}ec}{B}\frac{1}{1 + \nu_{in}^{2}/ \varOmega_{i}^{2}} \\ \varDelta = \bigl( 1+3\cos^{2}\theta \bigr) ^{1/2} \qquad b_{s}=\bigl(r_{E}^{3}/r^{3}\bigr) \varDelta \qquad s=qr_{E} \end{gathered}$$

Here, $g_{p}$ is the gravity in the $p$ direction, $\varOmega_{\alpha }= eB/m_{\alpha }c$ is the cyclotron frequency of species $\alpha$, and $\nu_{\alpha n}$ is the neutral collision frequency with the $\alpha$ species.