Self-consistent Global Transport of Metallic Ions with WACCM-X

The NCAR Whole Atmosphere Community Climate Model with thermosphere and ionosphere eXtension (WACCMX) v2.1 has been developed to include the neutral and ion-molecule chemistry and dynamics of three metals (Mg, Na, and Fe), which are injected into the upper mesosphere/lower thermosphere by meteoric ablation. Here we focus on the self-consistent electrodynamical transport of metallic ions in both the E and F regions. The model with full ion transport significantly improves the simulation of global distribution and seasonal variations of Mg. Near the magnetic equator, the diurnal variation in upward 5 and downward transport of Mg is generally consistent with the “ionosphere fountain effect”. The thermospheric distribution of Fe is shown to be closely coupled to the transport of Fe. The effect of ion mass on ion transport is also examined: the lighter ions (Mg and Na) are transported above 150 km more easily than the heavy Fe. We also examine the impact of the transport of major molecular ions, NO and O2 , on the distribution of metallic ions.


10
The presence of layers of meteor-ablated metal atoms between 80 and 105 km has been known for decades .
More recently, there have been a growing number of observations of the thermospheric metal layers up to 200 km. For example: Chu et al. (2011) reported neutral Fe layer observations up to 155 km at McMurdo, Antarctica (77.8 • S, 166.7 • E); Gao et al. (2015) found that several observations of Na layers reached up to 170 km, using a large-aperture astronomical telescope at Lijiang,China (26.7 • N, 100.0 • E); and Friedman et al. (2013) investigated a descending thermospheric K layer up to ∼155 15 km at Arecibo, Puerto Rico (18.35 • N,66.75 • W). These high altitude thermospheric neutral metal layers are challenging to explain, since ablation occurs predominantly below these altitudes, which raises interesting questions regarding their formation mechanisms. Neutral metal atoms and their corresponding atomic ions are tightly coupled through ionization (via photoionization or charge transfer with ambient ions e.g. NO + and O + 2 ), and neutralization (via dielectronic recombination with electrons, or dissociative electron recombination if they have formed a molecular ion) . Thus, the vertical and horizontal transport of metallic ions is a key process for controlling the behavior of neutral atoms in the thermosphere (E and F regions). Interestingly, the reported occurrence of thermospheric metal layers appears to show great geographical variability, underlining the importance of global ion transport. In addition, metallic ions play a central role in the formation of thin, concentrated layers of ions in the E-region (sporadic E layers, or Es), which affect radio transmission (e.g., Narcisi, 1968;Layzer, 1972;Plane et al., 2015;Yu et al., 2021). 25 Thus far, a number of modeling studies have attempted to simulate the transport of metal ions in the thermosphere and, more recently, the role of metallic ions in the formation of the thermospheric metal layers. For example, Carter and Forbes (1999) developed a two-dimensional (2-D) model to examine both global and local transport of Fe + ions; Chu and Yu (2017)  in addition to a detailed description of the neutral and ion-molecule chemistry of these metals and the meteoric ablation source required to model the metal atom layers around 90 km. This is described in Section 2. Section 3 presents the findings of the model simulation, focusing on the seasonal and the diurnal variation of ions and the affects of ion electro-dynamical transport.
The final section includes a brief summary and a discussion of future directions with the model. The key chemistry and dynamical features are based on CAM4 and WACCM4 and are described in detail in Marsh et al. 50 (2013b), and Neale et al. (2013). Validated metal chemistry modules for magnesium (Langowski et al., 2015), sodium (Marsh et al., 2013a), and iron (Feng et al., 2013) with updated rate coefficients from Plane et al. (2015), Bones et al. (2016) andViehl et al. (2016), are added. The transport of the neutral and ionized metallic species by eddy/molecular diffusion and winds is treated in the same way as most active chemical species (for example, O 3 ,CO 2 etc.).
A detailed description of WACCM-X 2.0 is provided by Liu et al. (2018a) and a brief summary is given here. The model 55 has some key features and improvements since WACCM-X 1.0 (Liu et al., 2010), including a self-consistent electrodynamics module, F-region O + transport, a solver for electron and ion temperatures, and reduction in the damping coefficient of atmospheric tides (Liu et al., 2018b). The model top is set at 4.1×10 −10 hPa (∼500 to ∼700 km, depending on solar activity), with a vertical resolution of a quarter of a scale-height in the mesosphere and thermosphere. The horizontal resolution is 1.9 • in latitude and 2.5 • in longitude, and all model results used in the paper are from one-year free-running simulations perpetual 60 year 2000AD under solar medium conditions (constant F107=124 and Kp=2.17) with an output frequency of 1 hr. The O + transport method is described in detail by Liu et al. (2018a).

Ion Transport Equation
Since meteoric ablation, deposition, transport by the neutral winds, eddy and molecular diffusion, chemical production, and loss are already contained in the metal chemistry modules, the metal ion transport is calculated separately, in a similar way to 65 the treatment of O + transport described by Liu et al. (2018a). The continuity equation of metal ion transport can be simplified as: where n i represents the number density of metal ions, and V i is the ion transport velocity. The ion transport velocity is adapted from the derivation described by Carter and Forbes (1999) and Chu and Yu (2017), extended to a 3-D global model: where E and B are the electric field and the Earth's magnetic field, respectively. V n is the neutral wind, and ξ = νin ωi is the ratio of ion-neutral collision frequency (in the laboratory frame-of-reference (Banks and Kockarts (1973)) to the ion gyro-frequency.
V ambi is the ion velocity due to ambipolar diffusion, given by Schunk and Nagy (2000, equations (5.54) and (5.70)), where 75 we treat the metallic ions as minor ions and O + as a major ion species in the F region. The first two terms are related to contributions by the neutral winds, with the first term being the V×B drift i.e. the Lorenz force. The third and fourth terms are due to the electric field. A Flux-Corrected Transport (FCT) algorithm (Boris et al., 1993) is applied to compute the vertical transport velocity; this algorithm is designed for solving steep density gradients.
As mentioned above, charge transfer from molecular ions such as NO + and O + 2 to metal atoms provides a major sources of 80 metallic ions. However, due to the short lifetime of these two molecular ions (∼5 mins in the daytime), they were assumed to be in chemical equilibrium in WACCM-X 2.0. In this study, the transport of molecular ions NO + and O + 2 is now considered along with the metal ions. Simulation results with and without the molecular ion transport are compared to determine the impact on the metal ion distribution. Na + , Fe + and Mg + undergo similar transport forces in the E and F regions, apart from the effect of their mass differences.
Because Mg + is the only one of these ions for which near-global observations are available (Langowski et al., 2015), we focus here on the results for Mg + . However, we also explore Fe/Fe + ion-neutral coupling in the formation of thermospheric metal layers, since Fe is the most abundance and widely studied thermospheric metal species. In addition, the Mg + /Fe + and Mg + /Na + ratios are employed to examine the effect of metal ion mass difference, because Fe + is more than twice as heavy as 90 Mg + (56 versus 24 amu), while the masses of Na + and Mg + are very similar (23 versus 24 amu). Figure 1 shows the monthly mean density of Mg + as a function of latitude and altitude, where the monthly mean data is zonally averaged. An obvious seasonal signal is exhibited with clear latitudinal dependence, which is generally consistent with the SCIAMACHY measurements ( Figure 6 in Langowski et al. (2015)). At middle latitudes (∼ 40 • ±10 • ), the peak altitude 95 of Mg + is ∼10 km higher in the summer hemisphere, which is very similar to the observations (Langowski et al., 2015). This appears to be caused by the field-aligned transport of the ions driven by the summer to winter neutral wind (second term in Eq.

Seasonal variation of Mg + simulated by WACCM-X
2), with monthly mean upward drift velocity of ∼5 m/s at a height of 150 km in the summer hemisphere (contour line in Figure   1; the seasonal variation of the neutral meridional wind is not shown). Since the metallic ions are the main reservoir for neutral metal atoms in the lower thermosphere , this is in good agreement with the summer maximum occurrence 100 trends of thermospheric neutral sodium layers observed by mid-latitude lidars (Wang et al., 2012;Dou et al., 2013;Xun et al., 2020). In contrast to the SCIAMACHY measurements, which shows a minimum at the equator, the WACCM-X simulation shows a maximum in peak altitude and number density at the equator. Note that the SCIAMACHY observations are made at a particular local time of around 10:00 LT, whereas the WACCM-X data in Figure 1 is a diurnal and zonal average. To address this, we also present the simulation results at the same local time (10:00 LT) in Figure 2, and they are in better agreement 105 with SCIAMACHY observations. Another noteworthy feature is that this simulation shows the pronounced maximum in peak altitude and density at ∼45 • (N/S) at 10 LT, in accord with the SCIAMACHY observations, which is absent in the previous models.  Langowski et al. (2015)), and the SBUV nadir measurements (Figure 9 in Joiner and Aikin (1996)). The partial column density above 110 km still shows a similar seasonal variation, which is probably due to thermospheric ion transport (not shown).

120
To investigate the diurnal variation of metallic ions in the model, we present the Mg + number density (on a log scale) as a function of latitude and altitude at the December Solstice and 0 • longitude; the panels show universal times (i.e. local times) of 00, 06, 12 and 18 UT (Figure 4 (a-d)). The white dotted line denotes the F 2 -layer height of the peak electron density (hmF 2 ).
This shows that the strongest diurnal variations are found in equatorial and high latitudes. Here we focus on the "fountain effect" on ion transport, where the equatorial ions are first lofted to higher altitudes via E×B motion, and then drift down along 125 the magnetic field lines (Kelley, 2009). Mg + is expected to be lofted to high altitudes (∼ 400 km) by the E×B drift above the magnetic equator during the day, because of the daytime eastward electric fields (Huba et al., 2019).  Oct. (e)). The phase of Mg + diurnal variation shows a high correlation with variations in the electron density (the change of hmF 2 ).
Instead of being redistributed along the magnetic field lines to the subtropical region by the "classical" fountain effect (e.g., Pi et al., 2009), the Mg + shows a more complex downward trajectory, i.e., the ions are not transported symmetrically to both sides of the geomagnetic equator, which is closer to the scenario proposed by Cai et al. (2019).

Fe and Fe + vertical profile comparison 135
In order to demonstrate the effect of electro-dynamical transport on the metals, a standard simulation without metal ion transport was performed (termed the control run). Figure 5 shows the Fe + and Fe vertical profiles for the ion transport run (solid lines), and control run (dotted/dashed lines) at the equinoxes and solstices. Three geographic latitudes at a longitude of 180 • are chosen for comparison: 0 • for the magnetic equator, 20 • S for the subtropical region corresponding to the fountain effect downward drift, and 45 • S for the middle latitude corresponding to the summertime peak altitude. Without the ion transport (the control 140 run), both Fe and Fe + exhibit roughly Gaussian-shaped layers with peak heights between 90 and 100 km (dotted/dashed lines).
When ion transport is turned on, the vertical profiles of Fe + vary depending on latitude and season. For instance, Fe + near the equator is always transported to a higher altitude (blue lines), consistent with the fountain effect in the dip equator region.   However, ions at subtropical (20 • S) and middle (45 • S) latitudes exhibit quite different transport motions. At middle latitudes, at midday the distribution is similar to that in the control run, but at midnight Fe + is transported to a high level (∼1 cm −3 at 145 200-300 km) except at the June Solstice (green lines). In contrast, the ions at subtropical latitude (orange line) are transported upward to a small extent at midday, but transported downward to a lower height (∼1 cm −3 at ∼100 km) at March Equinox and June Solstice at midnight, which is in reasonable agreement with the downward drift of the fountain effect along the magnetic field lines at midnight.
Since the neutral atoms are not directly transported by the electromagnetic field, they are influenced by ions through the 150 recombination between ions and electrons. The last two panels in Figure 5 illustrate the Fe atom distributions. In general, the upward transport of Fe + does not significantly contribute to changes of Fe (at densities > 1 cm −3 ), so that the vertical distribution of Fe in the transport run is similar to that in the control run. However, there is an obvious increase in highaltitude Fe (i.e. above ∼140 km) in the equatorial region (blue line), corresponding to the upward transport of ions; and at the December solstice around midnight, the Fe number density above 150 km is much higher than that in the control run at all 155 southern latitudes. Fe + , Sept. Equinox, 12 LTFe + , Sept. Equinox, 00 LTFe, Sept. Equinox, 12 LT Fe, Sept. Equinox, 00 LT  Figure 6 compares the Mg + /Fe + ratio (left panel) and the Mg + /Na + ratio (right panel) as a function of height and latitude at the equinoxes and solstices. Note that changes in the ratios below 100 km are due to differences in the ion-molecule chemistries of the metals , which is not the focus here. There are several advantages in choosing these three metallic ions.

160
First, Fe + is more than twice as heavy as Mg + . Second, Fe + is much heavier, and Mg + is slightly lighter, than the mean mass of air molecules in the E region. Third, Na + has a comparable mass with Mg + . As expected, the lighter ions are transported above 150 km more easily than the heavy Fe + , so that the Mg + /Fe + ratio increases from ∼1 at 120 km to >2 above 150 km and to > 30 above 300 km. By the same token, the Mg + /Na + ratio shows very little change above 120 km.
The zonally-averaged Fe + /Mg + ratio below 200 km simulated by WACCM-X also accords with the limited available obser-165 vations (Dymond et al., 2003;Kumar and Hanson, 1980), which showed that the average Fe + /Mg + ratio is around 1.5:1. The present study also simulates the extreme variability of the Fe + /Mg + ratio above 300 km (as low as 1:50 in Kumar and Hanson (1980)). However, the unexpectedly large Fe + /Mg + ratio (∼10-50) reported by Dymond et al. (2003) is only captured at a few points about 150 km in our model (not shown). Interestingly, the striking differences between distinct TIFe and diffuse TINa reported by Chu et al. (2020) is thought to be related to mass separation. There is no question that more observations are 170 needed to confirm and validate these findings.  In general, there is a good correspondence between the two simulations in terms of the latitude-altitudinal distribution of the monthly mean Mg + density. As seen in Figure 7c, the peak density in CE simulation is generally a little higher than that in 180 the TA simulation, especially in the high latitudes of the southern hemisphere. Figure 8 compares the diurnal variation of Mg + in the two simulations. Both cases simulate the significant "fountain effect", which was discussed in Section 3.2. With the transport of major molecular ions, the peak height of the metal layer after midnight is higher. As discussed by Plane et al. (2015), charge transfer of neutral metal atoms with NO + and O + 2 is the major sources of metallic ions in the E region. Due to the very short lifetimes of NO + and O + 2 during daytime, the transport of these molecular ions between model grid-boxes 185 has little effect on the metallic ions. In contrast, the reduced densities of the molecular ions (and electrons) at night means that their increased lifetimes become comparable to transport lifetimes. Additional metallic ions are therefore produced via charge transfer with the downward transport of NO + and O + 2 at night in the TA simulation.

Conclusions
The WACCM-X high altitude chemistry-climate model has been developed to incorporate the full life cycle of multiple me-190 teoric metal ions and atoms (Mg, Na, and Fe, currently). A major advantage of WACCM-X is the self-consistent treatment of dynamics and electrodynamics allowing us to quantitatively investigate the global distribution of metal ions and the formation mechanisms of thermospheric metal layers. The present study explores, for the first time, the seasonal variations of thermospheric metal ions by including global metal ion transport in the E and F regions.
There are a number of interesting findings: (1) A clear seasonal cycle is found in the monthly averaged global distributions 195 of Mg + , in good agreement with the SCIAMACHY measurements (Langowski et al., 2015). (2)   summer maximum occurrence of thermospheric sodium layers observed by mid-latitude lidars (Wang et al., 2012;Dou et al., 2013;Xun et al., 2020).
(3) Upward transport of metallic ions by E×B forcing is generally consistent with the "fountain effect".
(4) The formation of thermospheric neutral metal layers is strongly influenced by the upward transport of ions, since metallic 200 atoms and ions are coupled by relatively fast reactions in the lower thermosphere . (5) A pronounced mass separation of Fe + with the two lighter ions, Mg + and Na + , is demonstrated above 150 km, with the ratio between the lighter ions (Mg + and Na + ) and heavier ions (Fe + ) increasing with height by more than a factor of 2 above 150 km. More satellite observations of the Mg + /Fe + ratio are needed to test this prediction. (6) The role of NO + and O + 2 transport in the distribution of metal ions in the model is examined by comparing the two simulation results. It is found that they have little effect on the 205 monthly means of metal ions but affect the peak heights of metallic ions in the descending phase of the "fountain effect".
Previous research has established that thermospheric neutral metal layers are modulated by dynamics (e.g. gravity waves, atmospheric tides) (e.g., Chu et al., 2011;Xue et al., 2013;Cai et al., 2017;Qiu et al., 2016;Chu et al., 2020). In the future, this new version of WACCM-X can be used to investigate the effect of lower atmospheric dynamical processes on the formation of thermospheric neutral metal layers, by using the "specified-dynamics" version of the model (SD-WACCM-X).