First simulations of day-to-day variability of mid-latitude sporadic E layer structures

We present the first simulations that successfully reproduce the day-to-day variability of the mid-latitude sporadic E (Es\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_s$$\end{document}) layers. Es\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_s$$\end{document} layers appearing in the lower ionosphere have been extensively investigated to monitor and forecast their effects on long-distance communication by radio waves. Although it is widely accepted that the atmospheric tides are important in generating the Es\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_s$$\end{document} layers, no simulations to date have reproduced the Es\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_s$$\end{document} layers observed on a certain day. This is due to the lack of the combination of realistic information on the atmospheric tides in the lower ionosphere and a three-dimensional numerical ionospheric model that can simulate the precise transport of metallic ions. We developed a numerical ionospheric model coupled with the neutral winds from the GAIA (Ground-to-topside model of Atmosphere and Ionosphere for Aeronomy). The fundamental structures and the day-to-day variations of the Es\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_s$$\end{document} layers observed by a Ca+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Ca}^+$$\end{document} lidar are well-reproduced in the simulations.


Introduction
A sporadic E ( E s ) layer is a thin layer appearing in the lower ionosphere that consists of metallic ions such as Fe + , Mg + and Ca + . E s layers have been studied extensively since the middle of the twentieth century because they cause anomalous propagation of radio waves utilized for communication and broadcast (e.g., Whitehead 1970). More recently, it became known that the E s layers have a potential impact on air-navigation applications using very high frequency waves (Sakai et al. 2019).
It is widely accepted that the vertical shears of the horizontal winds of the diurnal/semi-diurnal atmospheric tides, which are dominant in the lower ionosphere, are responsible for the formation of the E s layers (Whitehead 1961;Axford 1963). In consideration of this, the E s layers are also called the tidal ion layers (Mathews 1998). According to the meridional wind-shear mechanism, the wind shears characterized northward winds above and southward winds below (for the northern hemisphere) accumulate ions to wind-shear nulls through the ionneutral collisions and ion motions along the magnetic field lines. According to the zonal wind-shear mechanism, the wind shears characterized westward winds above and eastward winds below accumulate ions to wind-shear nulls through the ion-neutral collisions and the Lorentz force exerted by the northward geomagnetic field. It is generally considered that not the meridional wind but the zonal wind shears cause the formation of E s layers especially below 120 km altitude (e.g., Haldoupis 2011). An E s layer is formed at the null and descends with the downward phase propagation of the tides. The typical descent velocity of the layer formed by the semi-diurnal tide is 0.6-3 m/s (Whitehead 1989).
Several simulations have been performed to investigate the tidal effects in E s layers (Carter and Forbes 1999;Krall et al. 2020;Mathews and Bekeny 1979). However, the neutral winds used in these simulations did not consider the day-to-day variability of the tides and thus, the E s layers were neither quantitatively nor qualitatively reproduced. Recently, a self-consistent, whole atmosphere-ionosphere model GAIA (Ground-to-topside model of Atmosphere and Ionosphere for Aeronomy) has reproduced 3-D thermospheric neutral winds exhibiting the day-to-day variability of realistic atmospheric tides by nudging the meteorological reanalysis data at the lower atmosphere (Jin et al. 2012). Shinagawa et al. (2017) estimated the occurrence rate of E s layers from ion convergences derived from the horizontal neutral winds of GAIA and compared it with the occurrence rate of E s layers deduced from GPS-CHAMP radio occultation (Shinagawa et al. 2017). The estimated occurrence rate of E s layers in general agrees with the seasonal and geographical variation of the observed rate. However, the timealtitude structures of the E s layers and their day-to-day variability have not been investigated because the exact transport process of metal ions has not been included in GAIA.
In the present study, we developed a 3-D numerical local ionospheric model for the metal ion layer with the neutral winds of GAIA. We simulate the fundamental structures and the day-to-day variations of E s layers, and then compare the variations with lidar observations performed at Tokyo.

Model description
We developed a 3-D numerical regional model of the ionospheric E region and the lower F region to study the three-dimensional structure of the E s layers. Ca + is used in the model as a proxy for metal ions. This model focuses on metal ion transport because the metal ions are long-lived (Murad and Williams 2002) and their chemical reactions are relatively unimportant. The continuity and momentum equations of the ions are written as where i represents the ion species such as NO + , O + 2 , O + , or Ca + , and e represents electrons. N is the number density, V the velocity, P and L the chemical production and loss, respectively, q the elementary charge, B the geomagnetic field, g the gravitational acceleration, the Boltzmann constant, m the mass, T the temperature, ν in the ion collision frequency with neutrals, and U the velocity of neutrals. ζ is the ratio of the ion collision frequency with neutrals to the ion gyrofrequency. ν in is obtained from Schunk and Nagy (2009). The chemical reactions in the model are shown in Tables 1 and 2.  Table 1 presents the chemical and photo chemical reactions of non-metal ions, and Table 2 those of Ca + . N + and N + 2 are included only in the chemical processes. T n = T i = T e is assumed in the model. A loss term due to three-body reactions involving Ca + is included below  Ferguson (1973) Schunk and Nagy (2009) Schunk and Nagy (2009) Schunk and Nagy (2009) Schunk and Nagy (2009) Richards et al. (1994) Richards et al. (1994) N + hv → N + + e − EUVAC model Richards et al. (1994) Richards et al. (1994) 100 km altitude in the continuity equation. The term increases with decreasing altitude and is given by The influence of the neutral atmosphere on the formation of E s layers is investigated in this study. The external electric field is not included in this model because the formation of E s layers is mainly caused by neutral wind-shears (Whitehead 1961;Axford 1963).
The numerical model is based on spherical coordinates. The positive directions are taken as vertically upward, southward, and eastward. The altitudinal coverage is from 85 km to 220 km. The grid spacing is set to the highest resolution of 0.5 km between 85 km and 140 km. Above 140 km, it gradually increases from 0.5 km to 2 km. The longitudinal and latitudinal domain size is 20 • × 20 • with a grid spacing of 0.25 • , and the center of the domain is located at 35.7 • N and 139.8 • E where the lidar observations were performed. The boundary conditions are set to zero gradients for all variables; however, at the top boundary, the densities of the ions except for O + are assumed to be in diffusive equilibrium. The number densities of neutrals except Ca and the neutral temperature are constrained by the NRLMSISE-00 model (Picone et al. 2002), and the geomagnetic field is taken from the IGRF-12 model (Thébault et al. 2015). The Ap and F10.7 index data are obtained from the World Data Center for Geomagnetism in Kyoto and NOAA/NCEI, respectively. The initial values of Ca and Ca + in the simulations are set to the annual mean profiles obtained by Plane et al. (2018), and the profiles are uniformly distributed in the horizontal domain. In order to precisely simulate ion transport, the Constrained Interpolation Profile scheme is used for the advection term (Yabe et al. 1991). The time step for solving equation (1) is 1 s. Prior to simulations, the chemical processes were calculated for 4800 s to obtain the chemical equilibrium state of each ion species. Then, the presented simulations were run for 36 h and the first 12 h were discarded to remove initial transients.
In the simulations in this study, the neutral wind velocity was obtained from GAIA. The horizontal resolution of the neutral atmospheric part of GAIA is T42 (the maximum horizontal wave number is equal to 42), corresponding to a grid spacing of 2.8 • longitude by 2.8 • latitude, and the vertical resolution is 0.4 scale height above the tropopause (e.g., Fujiwara and Miyoshi 2010;Miyoshi and Fujiwara 2008). The neutral wind velocity of GAIA was linear-interpolated in the vertical direction and spline-interpolated in the horizontal direction when it was used as an input for the regional ionospheric model. The diurnal and the semi-diurnal tides, which were extracted from the averaged 3-D neutral winds of GAIA in June, 2015 through a Fourier decomposition, were used as inputs for Case A to reproduce the fundamental structures of the E s layers due to the atmospheric tides. The F10.7 and AP index are set to averaged value in June, 2015. For Case B, the 3-D neutral winds at days corresponding to the lidar observations were used as inputs for comparison with the observed temporal Ca + profile variation. The F10.7 and AP indices are set to the values on each observation day. The lidar observations were performed at the National Institute of Polar Research, Tachikawa (Ejiri et al. 2019a, b). The time and altitude resolution of the observations were reduced to 600 s and 0.5 km, respectively, to match the output interval and altitude resolution of the model. The lidar observations have been operated only on eight nights in 2014-2016 because the lidar was on a development stage and not able to perform continuous observations. After the development, the lidar has been brought to the Antarctic. We were able to simulate E s layers only on six out of the eight nights because GAIA data are lacking for two out of eight nights. Hence, the results on the six nights shown in the next section are all of the cases we simulated. Figure 1 presents the temporal variations for Case A in which the monthly averaged tidal winds were used as inputs. The vertical axis represents the altitude from 85 km to 150 km, and horizontal axis the local time for 24 hours beginning from morning. The simulated Ca + density is shown in Fig. 1a, and the zonal, meridional  Rutherford et al. (1972) NO + + Ca → Ca + + NO k 19 = 4.6 × 10 −15 Rutherford et al. (1972) Ca + + e − → Ca + hv k 20 = 3.8 × 10 −18 ( Te 200 ) −0.9 Shull and van Steenberg (1982) Ca + hv → Ca + + e − 2 k 21 = 5.0 × 10 −5 Plane et al. (2018) and vertical component of the neutral winds utilized as the inputs are shown in Fig. 1b-d, respectively. Positive values in Fig. 1b-d represent eastward, southward and upward, respectively. Two E s layers are seen in Fig. 1a; "Layer 1" occurs in the daytime and "Layer 2" occurs in the nighttime. The E s layers appear at zonal wind-shear nulls rather than meridional-wind or vertical-wind ones between 115 and 120 km at 8 and 20 LT. As shown in Fig. 1b, the semi-diurnal tide is dominant between 115 and 120 km where Layer 1 and Layer 2 appear, and the descending layers are produced mainly due to zonal wind shears of the semi-diurnal tide. The descending velocity of the layers becomes smaller than that of tidal phase propagation due to the increase in the ratio of the ion collision frequency with neutrals to the ion gyrofrequency as the layers descend. Therefore, the E s layers gradually remain behind the wind shears. Subsequently, the E s layers stay around 100 km, and the descending layers merge with the stagnated layer. As a result, most of the Ca + accumulates around 100 km as shown in Fig. 1a. This numerical result is similar to previous E s layer simulation results (Carter and Forbes 1999;Krall et al. 2020;Mathews and Bekeny 1979). In Case A, the fundamental structure of the E s layer is reproduced; the E s layers are created by the semidiurnal tide, descend with its vertical wind shears, and remain around 100 km. In the observations, however, the E s layers have more complicated structures and show large day-to-day variabilities (e.g., Ejiri et al. 2019a, b;Raizada et al. 2020). The daily neutral winds of GAIA were used as inputs in Case B to reproduce and investigate the day-to-day variability of the E s layer structures. The simulated temporal variations of the Ca + altitude profiles are compared with those observed by a resonant scattering lidar on six nights from 2014 to 2016. Figure 2   In general, the simulation successfully reproduced the observations. In Fig. 2a and b, Layer 2 appears at higher altitudes and decays away after 4 LT. Layer 1 remains at the same altitude and disappears around 0 LT. On the other hand, two E s layers appear at 100 km and 110 km in Fig. 2c and d. The lower layer descends slightly and cannot be clearly seen after 4 LT. The altitudes at which Layer 1 and Layer 2 appear in Fig. 2e and f are similar to those of the layers shown in Fig. 2c and d, but the layers have no clear boundary between them and merge with each other at around 100 km. Only a single layer is created by the tides in each of Fig. 2g-l. The lifetimes of the layers also seem to be well-reproduced in Fig. 2g-l. As shown in Fig. 2f, h, and j, additional layers appear between 110 and 120 km in the morning. These layers do not appear in the simulations for the corresponding days, but a similar layer is reproduced at 110 km and 4-6 LT in Fig. 2a. The additional layer in Fig. 2l, which appears above 110 km and descends at a faster velocity than the downward phase propagation of the semi-diurnal tides, is also not reproduced in the simulation.

Discussion and conclusions
We have developed a regional ionospheric model of the metal ion layer to reproduce the day-to-day variability of E s layers. The regional model reproduced the day-to-day variation of the E s layers observed by the Ca + lidar well. In Case A in which the monthly averaged tidal winds were used as inputs, E s layers appear at ∼ 115 km altitude and descend to lower altitudes twice a day. The layers finally remain at ∼ 100 km altitude because the ions collide with neutrals more frequently below 100 km, resulting in the ion vertical velocity becoming almost equal to the slow vertical neutral wind velocity ( ≤ 0.1 m/s) in Case A. In contrast, in Case B where the daily neutral winds from GAIA were used as inputs, the E s layers on each night have different numbers, altitudes of appearance, and layer decay from the structures in Case A. In most of the nights shown in Fig. 2, Layer 2 appears and persists through the night. However, Layer 2 does not appear and Layer 1 persists till the morning in one out of the six nights as shown in Fig. 2g and h. The altitudes at which E s layers appear in Case B vary from 110 km to 130 km in the simulations and from 85 km to 140 km in the observations. It was found that the number of E s layers and the altitudes that they appear at depend on the strength and altitude of vertical zonal wind shears rather than the distribution of Ca + in the simulations. The dayto-day variability of the tides, which produces the variability of the vertical wind shears, is believed to be due to interaction between the tides and the planetary waves (Forbes et al. 2008;Liu et al. 2014;Miyoshi and Fujiwara 2003). The daily variation of the E s layer structures can be regarded as a manifestation of the day-to-day variability of the tides.
The E s layers that descend to ≤ 100 km disappear mainly because of the three-body reaction of Ca + . Whether or not the E s layers can descend to ≤ 100 km depends on the strength of the vertical neutral winds around 100 km altitude. Figure 3 shows the temporal variations of (a) vertical ion velocity, (b) zonal winds, (c) meridional winds and (d) vertical winds on December 18, 2015, where the E s layer descends clearly to < 100 km. The horizontal axis is same as Fig. 2c, but the range of altitude is from 85 km to 130 km altitude because the E s layers in Fig. 2c appear below 130 km and the vertical ion velocity above 130 km is mostly controlled by the meridional winds that is same as the previous numerical results (Carter and Forbes 1999;Krall et al. 2020). Positive values in Fig. 3a-d represent upward, eastward, southward and upward, respectively. The temporal variation of vertical ion velocity at 100-130 km has similar feature to that of zonal winds. In contrast, the feature of the temporal variation of vertical ion velocity below 100 km resembles that of vertical winds and not meridional winds. The vertical ion velocity below 100 km, which is equal to the vertical neutral wind velocity there, reaches up to ∼ 0.5 m/s in Case B and is larger than that in Case A. When the vertical neutral wind velocity is nearly 0.5 m/s, the layers in Case B descend to ≤ 100 km altitude and become neutralized by the three-body reaction of Ca + . The simulated Layer 1 shown in Fig. 2a, c, and e appears at a slightly higher altitude than the observed Layer 1. This indicates that the strength of the actual downward neutral winds might be stronger than that in the simulations. Previous numerical studies have neglected the influence of vertical neutral winds on E s layers (Carter and Forbes 1999;Krall et al. 2020;Mathews and Bekeny 1979), but many E s layers have been observed below 100 km (Ejiri et al. 2019a, b;Raizada et al. 2011Raizada et al. , 2020. Krall et al. (2020) have suggested that meridional winds account for the E s layers below 110 km because the SAMI3 simulations in "meridional winds only" case have reproduced observed layers below 110 km and the simulated layers have been associated with wind shears of meridional winds. As shown in Fig. 3, however, vertical ion motions below 110 km mostly result from zonal or vertical winds not from meridional winds. The E s layers below 110 km are not completely restricted by the wind shears (e.g., Haldoupis 2012). Therefore, E s layer structures might be influenced by horizontal transport by horizontal neutral winds rather than ion vertical convergence by the winds. Thus, it is crucial that the vertical neutral wind effects on E s layer structures below 100 km are included in simulations and the 3-D dynamics of simulated layers are investigated.
In addition to Layer 1 and Layer 2, other layers are seen in Fig. 2. They are morning E s layers that appear around 110-120 km and 3-5 LT as shown in Fig. 2f, h and j, and an E s layer that descends between 120 and 140 km with greater velocity than the phase velocity of the semi-diurnal tides shown in Fig. 2l. These discrepancies are induced because the neutral winds of GAIA do not contain all the atmospheric waves that cause E s layers because of the low spatial resolution. Figure 2a shows that the morning E s layer is reproduced in the simulation, but it is not observed as shown in Fig. 2b. In contrast, the observed morning E s layers in Fig. 2f, h and j are not reproduced. The simulated morning layer is a part of the 3-D structure of E s layers that moves horizontally. The 3-D structures of E s layers and their dynamics may play important roles in the complexity of the observed E s layers such as the morning layers and the fast descending layer. Further simulations are necessary to elucidate the effects of the atmospheric tides and, especially, the atmospheric gravity waves on the 3-D E s structures in terms of ionosphere-thermosphere coupling. In the present study, we emphasize the effect of the neutral winds on the E s layers because it is known for the neutral winds to affect the fundamental structures of the E s layers (e.g., Whitehead 1989;Mathews 1998;Haldoupis 2012), and the effect of the polarization electric fields on the E s layers is beyond Fig. 3 Temporal variations of (a) vertical ion velocity, (b) zonal winds, (c) meridional winds and (d) vertical winds on December 18, 2015. The horizontal axis is same as Fig. 2c-d, but the range of altitude is from 85 km to 130 km altitude and differs from that in Fig. 2