Numerical Simulations of Metallic Ion Density Perturbations in Sporadic E Layers Caused by Gravity Waves

Tidal signatures in sporadic E (Es) layer have been confirmed by observations and simulations. However, the effect of gravity waves (GWs) on the Es layer formation process has not yet been fully understood. In this paper, the modulation of Es layers by GWs is examined through numerical simulations, in which a physics‐based model of Es layer is forced by neutral winds from the High Altitude Mechanistic General Circulation Model that can resolve GWs with horizontal wavelengths longer than 156 km (λh > 156 km). Comparison of the simulation results with and without the GWs (1,350 km > λh > 156 km) forcing reveals that the inclusion of GWs leads to short‐period (1.2–3 hr) density perturbations in Es layers, which are also seen in ground‐based ionosonde observations. At a given time, the metallic ion density at altitudes between 120 and 150 km can either increase (by up to ∼+600%) or reduce (by up to −90%) in response to GW forcing. The relative density perturbations are smaller (by up to 60%) between 90 and 120 km altitude. It is also found that the GW effect on the metallic ion density relates to the longitude, which is mostly explained by the geographical distribution of GWs activity in the mesosphere and lower thermosphere region. The longitudinal variation of the background geomagnetic field plays only a secondary role.

• The modulation of Es layers by gravity waves (GWs) (1,350 km > λ h > 156 km) is simulated • GWs cause short-period (1.2-3 hr) Fe + density perturbations, leading to production of fine structure in mid-latitude Es layers • Longitudinal dependence in Fe + density perturbations is mainly explained by the longitudinal distribution of GWs activity in the mesosphere and lower thermosphere region

Supporting Information:
Supporting Information may be found in the online version of this article.
Correspondence to: T. Yu,yutaommnn@163.com Citation: Qiu, L., Yamazaki, Y., Yu, T., Becker, E., Miyoshi, Y., Qi, Y., et al. (2023). Numerical simulations of metallic ion density perturbations in sporadic E layers caused by gravity waves. Earth and Space Science, 10, e2023EA003030. https://doi. org/10.1029/2023EA003030 Huang and Kelley (1996) used a computer simulation to study GW modulation of existing Es layers and demonstrated that a horizontally stratified Es layer can be deformed by GWs and become a quasi-periodic wavelike structure. Using the non-dissipative linear GW model, Didebulidze et al. (2020) simulated the formation of multilayered Es and suggested theoretically that GWs can cause heavy metallic ion redistribution.
In the present study we use a 1-D Es layer model (Zuo et al., 2006) driven by neutral winds from the Ground-to-topside model of Atmosphere and Ionosphere for Aeronomy model (GAIA; Miyoshi & Fujiwara, 2003;Kobayashi et al., 2015), Whole Atmosphere Community Climate Model with thermosphere and ionosphere eXtension model (WACCM-X 2.1; H.-L. Liu et al., 2018), and High Altitude Mechanistic general Circulation Model (HIAMCM; Becker & Vadas, 2020;Becker et al., 2022) to examine the physical process of Es layer evolution, respectively. First, we compare the HIAMCM's capability to reproduce the tidal variation in Es layers with GAIA and WACCM-X that were used in previous studies (e.g., Andoh et al., 2022;Wu et al., 2021). Second, we focus on the physical process of Es layer evolution modulated by GWs, which is not taken into account in previous work.

Models and Observations
This work used an Es layer model (Zuo et al., 2006) and six ionosondes located at middle latitudes to simulate and observe the Es layer evolution process during the summer of 10-11 July 2007. The Es model solves for the altitude profile of the Fe + density as a function of time, based on the ion continuity and velocity equations. The ion continuity equation is given by where N i is the number density of Fe + , z is altitude, and q i and L i represent the chemical production and loss rates, respectively. Our simulations focus on the dynamical processes of Fe + , assuming chemical steady state (q i − L i N i = 0) (Didebulidze et al., 2020). The vertical velocity of Fe + , V iz is given by where γ is the ratio of the ion-neutral collision frequency (ν i ) to the ion gyrofrequency (ω i ) (Nygren et al., 1984). U x , U y , and U z are the southward, eastward and upward winds, respectively. The geomagnetic inclination θ is obtained from the International Geomagnetic Reference Field (IGRF11) (Finlay et al., 2010). The neutral winds are derived from atmospheric models, which will be described later. The initial density profile of Fe + is derived from the WACCM-Fe model, which is a global atmospheric model of meteoric iron (Feng et al., 2013). The ion continuity Equation 1 is solved by the finite difference method with a time step of 5 s and an altitude resolution of 500 m. The altitude range is from 80 to 150 km. The boundary conditions are set as ∂N i /∂t = 0.
The GAIA model is an upper extension of the Kyushu University General Circulation Model (Miyahara et al., 1993) that extends to the exobase at an altitude of ∼500 km. The neutral atmosphere part of GAIA assimilates meteorological reanalysis data below 30 km (Japanese 55-year Reanalysis; Kobayashi et al., 2015) by a nudging method (Jin et al., 2008). It has a grid of 2.8° longitude by 2.8° latitude horizontally. The GAIA model can simulate the large-scale waves with λ h > 1,000 km. The effect of GWs with λ h < 1,000 km on the background winds is taken into account by GW-drag parameterization (Miyoshi & Fujiwara, 2003;Sato & Yasui, 2018).
The WACCM-X 2.1 is a whole atmosphere model developed by the National Center for Atmospheric Research. It can self-consistently resolves the dynamic and physical processes from the surface to ∼700 km altitude (H.-L. Liu et al., 2018), which is also constrained by the MERRA-2 reanalysis data. The horizontal resolution is 1.9° in

Results
Figure 1 compares the simulated Fe + densities derived from the Es layer model driven by HIAMCM with total waves, HIAMCM with only large-scale waves, GAIA, and WACCM-X in panels (a-d), to confirm the HIAM-CM's capability to reproduce the tidal variation in Es layers and highlight what was not taken into account in previous work. The semidiurnal variation in Fe + density is reproduced through the Es layer model driven by HIAMCM, which is consistent with previous results (e.g., Andoh et al., 2022), suggesting that the Es layer (Fe + density >∼3.5 × 10 3 cm −3 ) occurrences are mainly controlled by neutral wind shears associated with atmospheric tides ( Figure S1 in Supporting Information S1). Compared with the simulations in Figures 1b-1d, Figure 1a shows some fine structures of Es layers, which is not simulated in previous results (e.g., Andoh et al., 2022; Wu . This is because HIAMCM can simulate GWs down to horizontal wavelengths of 156 km (Becker & Vadas, 2020), which drive the short-period perturbations in Fe + density within the Es layers. Note that the difference between Figures 1b-1d may be caused by differences in parameterizations of the physical processes, such as GW drag parameterization and solar radiation parameterization, and meteorological reanalysis data used in the atmospheric models. Figure 2 shows the temporal variations of simulated Fe + density and observed foEs over six mid-latitude stations during 10-11 July 2007. The Es layer model is driven by neutral winds from HIAMCM with total waves. The black solid lines represent the foEs recorded by ionosondes. First, both simulations and observations show that the Es layers tend to exhibit semidiurnal variations. Second, the longitudinal variations of simulated Fe + density also agree reasonably well with the observations. Third, short-period perturbations (1.2-3 hr) in the Es layers are displayed in both simulations and observations. Figure 3 shows the relative Fe + density perturbation (ΔFe + ) due to GWs (1,350 km > λ h > 156 km), over the six mid-latitude stations during 10-11 July 2007. The ΔFe + is given by (Fe + total waves − Fe + large−scale waves )∕(Fe + large−scale waves ) where the Fe + total waves is the simulated Fe + density with total waves, and Fe + large−scale waves is the simulated Fe + density with only the large-scale waves. On the one hand, GWs modulate the Fe + density more effectively at higher altitudes (above 120 km) than at lower altitudes. The Fe + density can be increased by 200%-600% or completely dispersed above 120 km, but changes by only about 60% below 120 km. On the other hand, Fe + density perturbations also show clear longitudinal dependence, being stronger over An Yang (127°E, 37°N) than other stations.

Discussion
The simulation experiments in this study suggest that GWs can modulate the intensity and structure of Es layers. The Fe + density shows short-period perturbations (1.2 hr at ∼135 km altitude, 3 hr at ∼100 km altitude) ( Figure S2 in Supporting Information S1), when GWs (1,350 km > λ h > 156 km) are included in the simulation. This means Figure 2. The temporal variations of simulated Fe + density and observed foEs (Es critical frequency) over six mid-latitude stations during 10-11 July 2007. The Es layer model is driven by neutral winds from High Altitude Mechanistic General Circulation Model with total waves (horizontal wavelength larger than ∼156 km). The black solid lines represent the foEs recorded by ionosondes. The simulations were run for 72 hr and the first 24 hr were discarded. • and ▿ on the horizontal axis represent the local noon (12:00 LT) and local midnight (24:00 LT), respectively. that upward-propagating GWs can cause sharp density gradient within Es layers. Miyoshi and Fujiwara (2008) examined the characteristics of GWs in the MLT region and found that the dominant periods of GWs for zonal wave number 80 (λ h ≈ 500 km) is 1.5-2 hr at 100 km altitude, which is basically consistent with the period of the Fe + density perturbations in our simulation results. The observations of foEs (Figure 2) also show perturbations at this period range.
The altitude dependences can be seen in the amplitude of ΔFe + . The ΔFe + is largest at 120-150 km (Figure 3). At lower altitudes (below 120 km), the GWs can still increase or reduce the Fe + density by about 60%. The modulation of GWs on metallic ion density is more effective at higher altitudes where the geomagnetically-driven vertical plasma transport is less hindered by ion-neutral collisions (Haldoupis, 2018). It is also worth mentioning that ion-molecule chemistry reactions (Plane et al., 2015) and ion ambipolar diffusion (Didebulidze et al., 2020) may affect the modulation of GWs on the Es layer intensity, which should be further examined through simulation and observation in the future. Besides, the Fe + perturbations moving downward with time, indicating an upward energy propagation of GWs. Above 120 km their phase velocity (or vertical wavelength) is very large, so they appear to have "quasi-vertical" structure. Figure 3 also shows that the amplitude of the GW-induced ΔFe + relates to the longitude. This is further highlighted in Figure 4a using the vertical ion convergence (VIC = −∂V iz /∂z, where V iz is the vertical velocity of Fe + ). It is known that VIC is that main physical parameter which controls the metallic ion density (Shinagawa et al., 2017). VIC is affected by both the wind shear caused by GWs and the geomagnetic field. Figure 4a shows the longitudinal distribution of the VIC perturbations due to GWs forcing (1,350 km > λ h > 156 km) at an altitude of 120 km and a latitude of 36°N. The variance of the VIC perturbations evaluated at each longitude is displayed in Figure 4b (blue line). The VIC variance shows a double-peak longitudinal structure. One local maximum is in the East Asian region (∼100°E) and another local maximum is in the American region (∼300°E).
The longitudinal dependence of the VIC variance can result from (a) the longitudinal variation of GWs activity and (b) the longitudinal variation of the geomagnetic field. In what follows, the relative importance of the two is examined. Figure 4c shows the global distribution of the kinetic energy of GWs (E k ) at 120 km altitude during 10-11 July 2007, which is given by = 1∕2 (Tsuda et al., 2000), where u′, v′, and w′ are the zonal, meridional, and vertical wind perturbations caused by GWs (1,350 km > λ h > 156 km). At mid-latitudes in the northern hemisphere, there are two GWs energy enhancements over East Asia and America, which corresponds to the two peak regions of the VIC variance. Miyoshi et al. (2014) investigated the global view of GWs activity (λ h > 380 km) in the lower thermosphere. At mid-latitudes in the northern hemisphere, there are two regions with enhanced GWs activity near ∼100°E and ∼300°E, which agrees with our results (Figure 4c).
The geomagnetic field configuration also affects the VIC (Haldoupis, 2011;Whitehead, 1960). In Figure 4b, the red line represents VIC variances as the blue line but the intensity B and inclination θ of the geomagnetic field are replaced with their corresponding zonal mean values (blue and red dashed lines in Figure 4d) so that the geomagnetic field is longitudinally invariant. The difference between the blue and red lines represents the effect of the inhomogeneous geomagnetic field. The difference is relatively small, suggesting that the longitudinal dependence of the VIC variance is primarily caused by the longitudinal variation of GWs activity rather than the longitudinal variation of the geomagnetic field. Nevertheless, the effect of the geomagnetic field on the VIC variance is visible, especially in the East Asian (∼120˚E) and American (∼270˚E) sectors. According to the zonal wind shear mechanism, ( ∕ )cos × ( 2 )∕( 2 2 + 2 2 ) , the stronger magnetic intensity B and lower magnetic inclination θ favor the formation of Es layer, where m i is the ion mass, e is the ion charge, and ν in is ion-neutral collision frequency. Figure 4d shows the longitudinal variations of B (blue solid line) and θ (red solid line) at a latitude of 36°N. In the East Asian sector, the strong magnetic field and low magnetic inclination contribute to an increase of the VIC variance, while in the American sector, the large magnetic inclination suppresses the VIC variance. The Es layer model in this paper successfully simulates the short-period variations within Es layers. However, the model does not simulate all the observed foEs peaks (Figure 2), firstly because variability of the electric field, turbulence, and sudden injection of meteoric ions all affect Es layer formation; and secondly, because our 1-D model does not capture the horizontal advection of metallic ions and inhomogeneity of the horizontal wind.

Conclusions
This work examined the physical process of Es layer evolution at middle latitudes during the summer of the 10-11 July 2007 by using an Es layer model driven by neutral winds from the HIAMCM model that includes the GWs (λ h > 156 km). The simulation results reproduce the tidal variations in Es layers, which agrees with the previous results from the Es layer model driven by GAIA and WACCM-X models. Besides, the inclusion of GWs in the simulations reveals fine structure in the Es layers caused by short-period (1.2-3 hr) GW-driven perturbations of the metallic ion density. The ion density can be increased by 200%-600% or completely dispersed between 120 and 150 km altitudes. The amplitude of ion density perturbations is about 60% between 90 and 120 km altitudes. In addition, metallic ion density perturbations also show longitudinal dependence, being stronger over East Asia (∼50-150°E), which is mainly explained by the geographical distribution of GW activity in the MLT region.

Data Availability Statement
The ionosonde data for this paper are available at: https://giro.uml.edu/didbase/scaled.php. The GAIA data are available at: https://gaia-web.nict.go.jp/data_e.html. The WACCM-X model source code can be found here (https://www2.hao.ucar.edu/modeling/waccm-x). The HIAMCM model simulations were performed by Erich Becker. Model documentations can be found in Becker and Vadas (2020) and Becker et al. (2022). The data used in the present simulations can been found in T. Yu (2023