Counterstreaming Cold H+, He+, O+, and N+ Outflows in the Plasmasphere

The Naval Research Laboratory (NRL) Sami3 is Also a Model of the Ionosphere (SAMI3) ionosphere/plasmasphere code is used to examine H+, He+, N+, and O+ thermal outflows during a storm. Here, H+ and He+ outflows are associated with refilling while O+ and N+ outflows are associated with ring current heating. An improved model of counterstreaming H+ outflows from the two hemispheres is presented, using an implementation of SAMI3 with two fluid species for H+. The two-fluid H+ model avoids nonphysical high-altitude “top-down refilling” density peaks seen in one-fluid H+ simulations. Counterstreaming cold ion populations are found in all cases. In these fully three-dimensional simulations with realistic magnetosphere boundary conditions, nonphysical top-down refilling density peaks were milder than those found in previous single-field-line or single-magnetic-longitude simulations. In the present two-fluid H+ case, “bottom-up refilling” density peaks were so mild as to be difficult to detect. For O+ and N+, the nonphysical high-altitude density peak is a brief (1–2 h) transient that occurs when heating-driven northward and southward flows first meet. In general, He+ outflows mimic H+ outflows while N+ outflows mimic O+ outflows.

One problem with fluid-code plasmasphere simulations is that the H + fluid tends to outflow from the north and south, colliding nonphysically near the apex (Banks et al., 1971;Richards et al., 1983;Singh et al., 1986). These outflows occur during the refilling of plasmasphere flux tubes (Sojka and Wrenn, 1985;Su et al., 2001;Dent et al., 2006;Sandel and Denton, 2007) following the erosion (Park, 1973;Foster et al., 2014) of the plasmasphere by a geomagnetic storm. In fact, the two outflowing proton streams normally pass through each other (Rasmussen and Schunk, 1988). In a single-fluid description of H + , the northward and southward velocities cancel when the two streams collide near the apex of the field line, producing a nonphysical density peak (Singh et al., 1986). While a comparison of fluid and semi-kinetic modeling (Singh et al., 1994) shows that this is a problem only during early-stage refilling, it is still a problem. Test runs (Krall and Huba, 2019) using the Naval Research Laboratory (NRL) Sami2 is Another a Model of the Ionosphere (SAMI2) ionosphere/plasmasphere code (Huba et al., 2000b) code show that a two-fluid description of H + dramatically improves modeling of earlystage refilling in comparison to the usual one-fluid description.
In this brief study, we revisit this issue using the NRL Sami3 is Also a Model of the Ionosphere (SAMI3) global ionosphere/ plasmasphere code (Huba and Krall, 2013). We consider SAMI3 simulations of a specific storm event, including plasmasphere erosion, refilling, and H + , He + , N + and O + outflows. In these simulations, O + (and N + ) outflows are driven by ring-current heating, as in Krall et al. (Krall et al., 2020). These results confirm past findings (Roberts et al., 1987;Craven et al., 1995) that the He + density generally mimics H + (Craven et al., 1997;Goldstein et al., 2003) and that N + generally mimics O + (Ilie and Liemohn, 2016).
We seek answers to the following two questions: does a twofluid description of H + improve SAMI3 modeling of early-stage refilling? and do single-fluid model O + outflows, like single-fluid model H + outflows, collide nonphysically?
We answer "yes" to both questions. In the simulations presented below, we show that two-fluid H + does indeed improve the modeling, but the effect is not as dramatic as that found in the test-case modeling of Krall and Huba 2019 and the prior single-field-line simulations of Banks et al. 1971 andSingh et al. (1986). We also show that, model O + outflows do collide nonphysically as outflows from the north and south meet each other at high altitudes. The result is a transient high-altitude density enhancement that fades within 1-2 h. Coincident N + outflows similarly collide, producing a brief transient. Figure 1 shows ground-based observations, Kp and Dst indices, for the October 7, 2015 (day 280) geomagnetic storm. Also shown is the extreme ultraviolet (EUV) F 10.7 index and its 81-days average, F 10.7A . These indices affect the strength of the ionosphere and the state of the background thermosphere. Figure 1 represents some of the inputs used in the SAMI3 modeling presented below. In addition to the EUV indices already described, the Kp index is used to determine the stormtime convection potential as in Volland (Volland, 1973), Stern (Stern, 1975) and Maynard and Chen (Maynard and Chen, 1975). The Dst index is used to determine ring-current heating of ionosphere and plasmasphere electrons as in Krall et al. (Krall et al., 2020).

SAMIRESULTS: HYDROGEN AND HELIUM IONS
We now present SAMI3 simulations of the 5-days pre-storm, storm, and recovery period shown in Figure 1. SAMI3 (Huba et al., 2005;Huba and Krall, 2013) simulates H + , He + , O + , N + , O + 2 , N + 2 , and NO + ions in the ionosphere and plasmasphere, with one fluid used for each ion. Below, we refer to this version of SAMI3 as either "SAMI3 (one fluid H + )" or, simply, "SAMI3." We also present simulations of this same time period using SAMI3 (two-fluid H + ), also called "twostream SAMI3." As in two-stream SAMI2 (Krall and Huba, 2019), two-stream SAMI3 represents the H + ion component of the ionosphere and plasmasphere using two model fluids, one with a source only south of the magnetic equator, the other only north of the magnetic equator. During quiet times, these two fluids combine diffusively to mimic the model H + density of SAMI3 (one fluid H + ).
The storm simulated here is the same as that in Krall et al. (2020). To ensure that the only differences between SAMI3 and two-stream SAMI3 reflect an updated treatment of the H + ions, we began by repeating that previous work with the latest version of the SAMI3 code. We then modified the code to implement two-fluid H + and repeated the simulation. As expected, the new SAMI3 run of the same storm is virtually identical to that presented in Krall et al. (2020). Because, as we shall see below, the effect of counterstreaming is not dramatic, results from the new two-stream SAMI3 run are also quite similar to those presented in Krall et al. (2020). However, our attention in that prior work was on O + , not H + .
Let us now turn our attention to H + and He + . Figure 2 compares SAMI3 (left column) to two-stream SAMI3 (right column) at identical times. Shown are color contours of electron density n e at longitude 0. Each plot is near midnight local time, where "empty" flux tubes convect in from the tailward boundary. Based on model-data comparisons at geosynchronous orbit (Krall et al., 2018), the H + density boundary value is set to  is caused by colliding H + streams from the two hemispheres. While a weaker version of feature A seems to be present in the two-stream SAMI3 result, panel (F), analysis shows that panel (F) also includes a second, much weaker peak. Two days later, panels (D) and (H), SAMI3 and two-stream SAMI3 come into closer agreement. 0.1 cm −3 . These times are chosen to include the peak of the storm, panels (b) and (f), where and when counterstreaming H + outflows from the ionosphere are strongest. Figure 2, panel (b), shows the one-fluid colliding upflow effect, with a prominent peak in the electron density near the apex of the geomagnetic field for altitudes > 3R E . This peak, labeled "A" in panel (b), appears to be significantly reduced (and shifted slightly southward) in the two-stream SAMI3 result, panel (f). This will be examined further below. Two days after the storm, panels (d, h), SAMI3 and two-stream SAMI3 give nearly identical results. Because these plots are near local midnight, the topside ionosphere electron hole described by Huba et al. (2000a) is evident at altitude 0.5 R E in each panel.
We now consider plots versus position along an L 5 field line at the same fixed longitude, where L is the McIlwain parameter (McIlwain, 1961). Figure 3 shows the colliding-upflow effect, with a peak in the H + density, panel (a), at latitude −10°. This is the same as feature A in Figure 2B. The corresponding velocity, panel (b), shows a region of near-zero H + velocity between high-speed (8 km/s) upflows from the north and south. Here, a positive velocity is northward.
Dashed lines in Figure 3 indicate He + . These show that, as expected from past modeling and observations, the He + density is generally coincident with the H + density. Interestingly, this coincidence includes the nonphysical, but likely insignificant, He + density peak at latitude −12°in panel (a). Similar to the way that the He + density mimics the H + density in panel (a), the He + velocity mimics the H + velocity in panel (b).
Two-stream SAMI3 results are shown in Figures 3C,D. Here, separate south and north components of H + are shown as dotted curves. In panels (c, d), we see that H + from the north (green dotted curve) moves southward at a high velocity. The velocity falls to zero as it moves southward. As a result, H + ions from the north accumulate in the south, creating a peak at latitude −20°. This is the density peak, seen in Figure 2F that appears to correspond to the much stronger peak of Figure 2B. In fact this peak in H + N is one of two H + peaks, the other being a peak in H + S at latitude 25°. In the latter case, H + ions from the south accumulate in the north. The H + velocity shown in Figure 3D, a density-weighted average of the two separate H + fluid velocities, shows that the two-fluid model of H + affects both densities and velocities.  are color contours of O + density at longitude 0. In both cases O + is represented by a single fluid; in two-stream SAMI3, only H + is represented by two fluids. Figure 4 shows a weak, local O + density peak at nearly the same latitude as the H + peak of Figure 2B. This seems to suggest that converging O + outflows from the north and south are colliding at this point, but this is misleading. In fact, these plots are near local midnight, where O + velocities are often downward. Unsurprisingly, a comparison between the left hand panels (SAMI3) and the right hand panels (two-stream SAMI3) of Figure 4 shows that the two-fluid treatment of H + has little effect on O + . For this reason, we presently focus our analysis of O + dynamics on the SAMI3 results. Figure 4C shows significant O + density at L > 3.5, with densities exceeding 100 cm −3 . This is the O + "shell/torus" (Horwitz et al., 1986;Roberts et al., 1987) that was first identified (Chappell, 1982) using the Retarding Ion Mass Spectrometer (RIMS) instrument on the Dynamics Explorer (DE) spacecraft. This aspect of these model results was discussed at length in Krall et al. (2020) [note: Krall et al. (2020) incorrectly states that the heating was applied only at high altitudes; in fact the heating was applied along the entirety of each field line; this issue will be explored further in the near future].

SAMI3 RESULTS: OXYGEN AND NITROGEN IONS
To discern the source of the high-altitude peak in the O + density, Figure 4B, we examine the dynamics along an L 4 field line at the time of Figure 4B and at a time 14.5 h later, when ring-current heating of ionosphere and plasmasphere electrons is driving significant O + outflows. Figure 5 shows O + and N + densities and velocities along an L 4 field line at these two times. In the afternoon sector and soon after the peak of the storm, Figures 5C,D, the ring current heating of plasmasphere electrons via Coulomb collisions with ring current ions is strongest. At this time the tendency of O + to be contained by gravity is overcome and upflows from the north and south are evident in panel (d). The converging velocity pattern and the local jump in the O + and N + densities at latitude 22°in Figures 5C,D suggest that these outflows might be better described using two fluids, as with H + in twostream SAMI3. This will be discussed further below.
After dusk, velocities in the south reverse, panel (b). Here, downward flows in the south are growing while upward flows in the north persist. Near the apex, the velocity is close to zero. Profiles of O + density and velocity in Figures 5A,B show that the density peak near the apex, visible only in panel (a), is not associated with simultaneous converging upflows. We speculate that the O + density peak near the apex of the field line in Figures 4B, 5A represents density that is slowest to descend at night. Dashed lines in Figure 5 show that N + generally mimics O + , but with a lower density.

DISCUSSION: BOTTOM-UP REFILLING
Based on past work, we expect a two-stream treatment of H + to replace "top-down" refilling, seen in one-fluid models (Banks   al., 1971), with "bottom-up" refilling (Rasmussen and Schunk, 1988). With top-down refilling, the initial peak is at high altitude. With bottom-up refilling, two initial peaks occur at lower altitude. In Figure 3, the top-down peak is in panel (a) at latitude −10°. In the two-stream SAMI3 result, panel (c), only one peak is clearly visible; it is located at −20°. In fact, there are two H + peaks in panel (c), the northern peak being at 25°. This can be seen in Figure 3, panel (c), where the separate H + S (H + coming from the south, blue) and H + N (green) densities are shown as dotted curves. These are the two peaks expected for bottom-up refilling.
In Figure 3C, H + S , has high values in the south (<−40°), where it is the main component of the total H + density, and in the north (25°), where it is contributing to the bottom-up refilling process. Similarly, H + N has high values in the north (>40°), where it is the main component of the total H + density, and in the south (−20°). The H + S and H + N field-aligned velocities, panel (d), confirm that H + S is rushing northward from the south while H + N is rushing southward from the north. In bottom-up refilling we expect the two peaks to form initially at lower altitudes (higher latitudes) and then move upward and equatorward as the flux tube fills (Rasmussen and Schunk, 1988) (Figure 3). Plots of H + S and H + N densities at other times, not shown, follow this pattern. It is notable that neither the one top-down peak, Figure 3A, nor the two bottom-up peaks, Figure 3C, are nearly as dramatic as the corresponding features found in test-case numerical models of refilling [e.g. 42, 23]. In fact, differences between these SAMI3 and two-stream SAMI3 results are confined to the midnight sector at high altitudes where the top-down density peak in SAMI3 exceeds the two-stream SAMI3 density by only a factor of 2. Otherwise, results from the two simulations are quite similar. While observations of low-energy ion outflows (Sagawa et al., 1987;Liemohn et al., 2005) and plumes (Borovsky et al., 2014;Lee et al., 2016) often show counterstreaming ion populations, corresponding "bottom-up" (or discredited "topdown") density enhancements have yet to be detected.
These runs illustrate the importance of multidimensional simulations for the purpose of accurately simulating refilling events. Specifically, our present runs differ from prior onedimensional (single field line) simulations [e.g. 42, 37] and our own prior two-dimensional SAMI2 (one magnetic longitude) simulations (Krall et al., 2008;Krall and Huba, 2019) in that the low-density pre-refilling conditions imposed in those lowdimensionality simulations are not present here. The "empty" field lines that are convected into a global plasmasphere simulation are not empty at all. They instead represent the cold plasma component of the tailward magnetosphere. In this case, we represented that plasma by imposing a minimum density of 0.1 cm −3 . Further, these inward-convecting field lines develop low-density counterstreaming ion flows as soon as they enter the simulation. While it might be possible to capture some of these features by imposing convection drifts and a tailward boundary condition on a SAMI2 simulation, a complete numerical description, including zonal convection drifts, are obtainable only via three-dimensional simulation. Ultimately, such modeling will need to be placed into the context of a global magnetosphere simulation, such as by Glocer et al. (2020).

DISCUSSION: COUNTERSTREAMING OUTFLOWS
To be clear, in this study of H + , He + , O + , and N + counterstreaming outflows, only H + counterstreaming is simulated, and only in the two-stream SAMI3 case. He + , O + , and N + are represented using a single fluid per ion species. In the latter cases, would-be counterstreaming flows are indicated by simultaneous converging upflows, such as in Figures 3B, 5D, with corresponding density peaks, Figures 3A, 5C.
From previous studies, we expect that the He + and H + components of the plasmasphere are generally co-located. For example, remote and in situ measurements show that the He + and H + components of the plasmapause are closely aligned (Goldstein et al., 2003). These results, specifically Figure 3, are consistent with this past understanding. They suggest that counterstreaming He + outflows, like H + outflows, occur during refilling (Richards et al., 1983). However, given the lack of dramatic differences between these SAMI3 and two-stream SAMI3 results, and that the He + density is typically less than 10% of the H + density, we speculate that counterstreaming He + outflows are not significant.
In the case of O + , we find simultaneous converging highaltitude upflows only during the day and during the storm, when the ring current heating outflow is further supported by the tendency of the ionosphere and plasmasphere to swell during the day (Galvan et al., 2008). In Figure 5D we find nearly identical converging upflow velocities for O + (solid curve) and N + (dashed curve). These converge at latitude 22°, where the field-aligned velocity crosses zero. At this point both the O + and N + densities jump by an order of magnitude and N + has a local density peak. The O + density curve at latitude 25°might also be interpreted as a local density peak. Figure 6 further examines the colliding outflow that is apparent in Figure 5D. O + density contours (left column) and corresponding O + and N + density profiles along the L 4 field line (right column) show that the density jump of Figure 5C, also shown in Figure 6F, occurs just as the north and south O + outflows meet at high altitude. Figures 6A,E, 30 min earlier, shows that these two relatively high density (100 cm −3 ) outflows have not yet connected. Figures  6C,G, 1 hour later, and (d, h), 2 hours later, show that the localized density peak dissipates very quickly. Velocity profiles (not shown) at these same times demonstrate that the coincidence of the velocity zero-crossing and the local density peak, see in Figures 5C,D, does not persist.
These results suggest that counterstreaming O + ions are both likely be present at high altitude during a storm and likely to be only a transient effect. We speculate that a two-fluid model of O + might not exhibit the same density jump as in Figure 5D; instead one might obtain a different, but similarly brief, transient density structure.
Unfortunately these O (1 eV) outflows are very difficult to measure. For example, Sagawa et al. (1987) uses Dynamics Explorer 1 data to simultaneously measure composition and pitch angle distribution, but only for energies >10 eV. Using these same data, Lin et al. (1982) infer counterstreaming Maxwellian ion populations such as simulated here. Lee et al. Frontiers in Astronomy and Space Sciences | www.frontiersin.org August 2021 | Volume 8 | Article 712611 8 (2016) uses measurements to infer numerous plume or outflow occurrences. However, because measurements used in that study cannot access ion energies below 5 eV, the cold plumes and outflows could not be described in detail. Nevertheless, even within the limitations of current observations, model/data comparison and analysis of well-measured events could provide a useful evaluation of these results.

CONCLUSION
As in past modeling studies (Richards et al., 1983), we find evidence for stormtime counterstreaming H + and He + outflows. In this global simulation of an actual storm, however, such counterstreaming is notable only in the midnight sector. Results also suggest the occurrence of stormtime, high-altitude O + and N + counterstreaming in the afternoon sector, where the ring-current heating that drives the outflow is strongest (Krall et al., 2020). In these simulations, where O + and N + are each modeled using a single fluid (only H + was modeled using two fluids), converging O + and N + produce small density peaks (less than an order of magnitude over background values), Figure 6F. These density enhancements dissipate rapidly, Figures 6G,H. In Figure 6F, the local N + density peak is on an L 4 field line at latitude 28°.
While two-stream SAMI3 (SAMI3 with a two-fluid model of H + ) leads to minor changes in the results at plasmasphere altitudes, the biggest effect is the avoidance of a small, nonphysical H + density peak near the equator, in the midnight sector, during early-stage refilling. At lower altitudes (<2R E ), this effect is less significant. For a strong storm, however, low-density flux tubes might convect in from the tail quickly enough to bring flux tubes to low L values while still undergoing early-stage refilling; this might produce stronger density peaks. In all cases simulated here, nonphysical high-altitude density peaks were quite mild, with the nonphysical density peak exceeding the expected result by less than a factor of 10. In this global simulation of an actual storm event, non-physical one-fluid-model density peaks were not nearly as pronounced as those seen in past refilling test-case studies (Singh et al., 1986;Rasmussen and Schunk, 1988;Krall et al., 2008;Krall and Huba, 2019).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10. 5281/zenodo.4771210.

AUTHOR CONTRIBUTIONS
Both authors, JK and JH contributed to this manuscript.

FUNDING
This research was supported by NRL base Funds, NASA Grand Challenge award NNH17AE97I, and NASA Living with a Star award 80NSSC19K0089.