Sporadic E morphology based on COSMIC radio occultation data and its relationship with wind shear theory

The S4max data retrieved from the Constellation Observing System for the Meteorology, Ionosphere, and Climate (COSMIC) radio occultation (RO) measurements during 2007 to 2015 is adopted to investigate the global distribution and seasonal variation of the sporadic E (Es) layers in the present work. The long-term and short-term global Es occurrence maps are presented and the spatial and temporal distributions of Es occurrence rates (ORs) are further confirmed and studied. The International Geomagnetic Reference Field model (IGRF12) is used to calculate the horizontal intensity and inclination of the Earth’s magnetic field. The analysis shows that the Earth’s magnetic field is one of the fundamental reasons for the global distribution of the Es layers. In addition, the Horizontal Wind Field model HWM14 and the IGRF12 model were employed to calculate the vertical ion convergence (VIC) to examine the role of neutral wind shear in the global distribution of the Es ORs. The results reveal that the middle latitude distribution of simulated vertical concentration of Fe+ is similar to that of Es ORs, which indicates that the VIC induced by the neutral wind shear is an important factor in determining the geographical distribution, summer maximum (or winter minimum) and diurnal characteristics of Es ORs in middle latitudes. The new findings mainly include the following two aspects: (1) in summer over mid-latitudes, VIC peaks in the morning and afternoon to evening, which explains the semidiurnal behavior of Es ORs; (2) VIC reaches its minimum value in low-altitude (100 km) areas, which is the reason for the significant decrease in Es ORs in low-altitude areas. The disagreements between the VIC and Es ORs indicate that other processes, such as the meteor influx rate, the ionospheric electric fields and atmospheric tides, should also be considered as they may have an important impact on the variation of Es layers.


Introduction
The ionosphere sporadic E (Es) layers are the ionized part at 90-130 km altitudes of Earth's upper atmosphere. Existing studies show that Es layers are thin-layer structures with a vertical thickness of less than a few kilometers and horizontal extension of up to several hundred kilometers with high densities of electron and metal ions (such as Fe + , Mg + , Ca + and Na + ) (Kopp 1977;Arras 2010). There is now enough evidence to suggest that sporadic E is not as "sporadic" as its name implies, but rather a regularly occurring phenomenon in middle latitudes (Haldoupis 2011), thus they may affect radio communications and navigation systems severely.
Es layers have been investigated extensively using the ground-based observations (Mathews 1998;Pignalberi et al. 2014;Resende et al. 2016). The ground-based radars have been used to observe the Es layers since decades ago. A global distribution map of Es layers was first presented based on ionosonde measurements (Leighton et al. 1962). The presence of tidal and planetary wave periodicities in Es layers were recognized based on critical frequencies measured by ionosonde stations in the middle latitudes (Haldoupis et al. 2004). Recently, the Global Positioning System (GPS) radio occultation (RO) measurements are widely employed to investigate Es layers (Hocke and Tsuda. 2001;Arras et al. 2008;Fytterer et al. 2014). Wu et al. (2005) used the GPS RO measurements from the Challenging Minisatellite Payload (CHAMP) mission to invert the global morphology of the Es layers. Arras et al. (2008) investigated the global distribution and seasonal variation of Es layers based on GPS RO measurements provided by the Constellation Observing System for the Meteorology, Ionosphere, and Climate (COSMIC) mission, and revealed the effect of the Southern Atlantic Anomaly (SAA) in Earth's magnetic field intensity on the reduction of the Es layer occurrence rates (ORs). Using the amplitude fluctuations of the GPS RO signal, Zeng and Sokolovskiy (2010) investigated the vertical extent of Es layers and showed that the thickness of Es layers ranged from 0.6 to 4 km with a most probable thickness of 1.5 km. Yu et al. (2021a) conducted the statistical comparisons between the Es layers derived from the COS-MIC RO scintillation index S4max and those derived from the critical frequencies (foEs) measured by ionosonde stations, and stated that the GNSS RO technique can provide global estimates of Es layers as a complement to ground-based observations.
Great efforts are made to explain the physical mechanism responsible for the formation of the Es layers (Whitehead 1989;Arras 2010;Haldoupis 2011). In the equatorial regions, the Es layers are well related to equatorial electro-jet, and the gradient-drift instability is the most popular theory to explain the formation of Es layers in these areas (Whitehead 1989;Tsunoda 2008). At polar latitudes, convection electric fields and the vertical motion of gravity waves are very efficient in concentrating the ions of Es (Kirkwood and Nilsson 2000;MacDougall et al. 2000a, b). In the middle latitudes, the wind shear theory is widely accepted to describe the physics mechanism of Es layer formation. It relies on the idea that thin layers of intense metal ions can form in the ionosphere region by vertical ion convergence (VIC) driven by vertical shears in the horizontal neutral wind, which was first proposed and formulated in the early 1960s by Whitehead (1961), andAxford (1963), and since developed further by other scholars (Nygren et al. 1984;Bishop and Earle 2003;Yu et al. 2019). Chu et al. (2014) employed COSMIC RO measurements and HWM07 model to study the relationship between the spatiotemporal distribution of Es ORs and that of the simulated divergence of vertical Fe + flux, and found that the seasonal variations of Es ORs are likely attributed to the convergence of the metal ion flux driven by vertical wind shear. Shinagawa et al. (2017) derived the VIC from the neutral wind simulations and suggested that the simulated VIC could partially explain the geographical and seasonal variations in Es ORs. Yu et al. (2019) compared the distributions of Es intensity and the divergence of vertical ion velocity, and found that wind shear can explain the seasonal dependence of Es intensity at the altitudes of 97-114 km, but it is hard to explain that at higher altitudes. Although plenty of work has been done about the Es morphology using GNSS RO data, mid-latitude regions are concerned about the most in previous studies, and detailed analyses about the distributions of RO-derived Es ORs over the polar regions and the related diurnal variations are still limited. Moreover, although the previous research results showed that the distribution and seasonal variation of Es layers can be roughly explained by the wind shears, the explanations of the diurnal variation and altitude distribution of the Es layers are still limited.
In summary, the overall distribution and the main seasonal characteristics of the Es layers have been basically explained, and its influencing factors include neutral wind shear, Earth's magnetic field, metal ion concentration, ionosphere, and atmospheric tides. However, the analysis and interpretation of the diurnal variation characteristics and the altitude distribution of Es layers have not been fully studied yet. Therefore, there are still a number of detailed issues related to the Es layer morphology which need to be resolved for a better understanding of the mechanism of Es layer formation. In the present study, we try to comprehensively analyze the relationship between Es ORs and their influencing factors, resolution of which may bring about a further progress in the study of the Es layers. Section "Data and processing" describes the data set of S4max indexes from the COSMIC RO mission and the corresponding data processing method. In the succeeding section, the spatial and temporal distributions of the Es ORs are presented. In the Sect. "Wind shear simulation", the spatial distribution and the time variation of the VIC are derived from the HWM14 model and the IGRF12 model based on the wind shear theory, then the relationship between the spatiotemporal distribution of Es ORs and the simulated VIC is discussed.
Finally, Sect. "Conclusions" points out the conclusions and the unresolved problems of the current research.

Data and processing
Scintillations in GPS signals occurring in the ionosphere are directly related to electron density fluctuations, and can be used to detect the ionospheric Es layers (Arras 2010;Chu et al. 2014;Yu et al. 2019). The COSMIC mission of the University Corporation for Atmospheric Research (UCAR) was launched in mid-April of 2006, and were still operating at the end of 2017 (Tsai et al. 2018). COSMIC consists of six low Earth orbit (LEO) satellites, which could observe more than 2000 RO events per day in the full operational stage of the mission, providing both 50 Hz and 1 Hz sampling rate observation data. The S4 index data are derived from the root mean square (RMS) signal-to-noise ratio (SNR) intensity fluctuation and mean SNR intensity in the received L1 GPS signal sampled at 50 Hz over the course of 1 s (Ko and Yeh 2010). Each RO event provides a vertical profile of the S4 index extending from Earth's surface to 800 km above Earth. The 1-Hz S4 scintillation index data are available from the 'scnLv1' files provided by the COSMIC Data Analysis and Archive Center (CDAAC). In the header of each 'scnLv1' file, the value of 'S4max' is denoted, and the parameters including the latitude, longitude, altitude, and time at the tangent point corresponding to the S4max value are also given. In the present paper, the study of the global distribution pattern of Es is carried out based on these S4max values. Figure 1 Fig. 1b, the number of available COSMIC RO S4max values in each year is greater than 1 million during 2007-2015. Since the amount of the COSMIC RO observation data decreases significantly due to the aging of the low Earth orbit (LEO) satellites during 2016 and 2017, the S4max data during the nine years from 2007 to 2015 are used to investigate the long-term global distribution pattern of the Es layers in the present study. Quality control on the raw S4max values is carried out considering that gross errors and outliers may exist in the raw data due to the intermittent anomalies in the C/A SNRs from antennas and observation noise (Rocken et al. 2000). According to the contents of the header of scnLv1 files, the valid range of the S4 index is between 0 and 5. Therefore, in the present study, after removing the null values in the raw S4max data set, we only select S4max data that fits the range of 0-5 as valid data. Approximately 20% of the raw S4max values are excluded by this screening procedure.
The altitude distribution of the valid S4max data is shown in Fig. 2a. It can be seen from Fig. 2a that the heights corresponding to the valid S4max values vary from 40 km to about 800 km, which is consistent with the previous study , and most of them are obtained at the altitude range between 40 and 120 km, with the peak frequency occurring at around 100 km and drastically low frequencies above 130 km. In addition, Fig. 2b shows the main spread of valid S4max values. In the subfigure, the red dotted line and the black line represents the average and the median of the valid S4max values, respectively, and the lower and upper lines of the rectangle represent the first quartile (25th percentile) and third quartile (75th percentile), whereas the lower and upper blue horizontal lines represent the 5th percentile and the 95th percentile of the valid S4max values. From Fig. 2b, it can be seen that about 90% of the S4max values are distributed between 0 and 1, and the percentage of the S4max values smaller than 0.4 is approximately 75%, which is consistent with the study of Yu et al. (2021a). Considering that the altitude range for the occurrence of Es layers is between 90 and 130 km, the second step involved in the preprocessing procedure is to extract S4max values within this altitude range, which is marked as qualified S4max values. In addition, the time variation of received signal power due to changes in transmission medium or paths is known as fading, which depends on various factors as path loss and multipath attenuation (Lee et al. 2017). If the RO signal experiences significant fading and scintillation at the same time, it has a great probability of having passed through the high-intensity ionosphere. Considering that signal intensity fadings do not occur for all scintillation events (Guo et al. 2019), it is mandatory to set a threshold for the S4 index. Only those S4 values greater than the threshold would be accompanied with obvious signal fading, indicating that the corresponding RO signals are more likely to have passed through the high-intensity ionosphere. Guo et al. (2019) analyzed the ratio of fading occurrence over scintillation events in relation to S4, and pointed out that the ratio increases from around 8% when S4 = 0.3 to nearly 80% when S4 = 0.5. This indicates a higher probability of signal intensity suffering from fading when S4 is greater than 0.5. Therefore, in the present work, a threshold of 0.5 is set for the S4max values, which means that the qualified S4max values between 90 and 130 km which are greater than 0.5 are designated as Es layer-related ones, and represent the intense Es layers. The threshold is also the same as those adopted by Ko and Yeh (2010) and Qiu et al. (2019).
The simplified process of data processing is shown in Fig. 3. After extracting all the qualified S4max values and the Es-related ones during 2007-2015, the qualified S4max values and the Es-related ones are gridded according to latitude-time, latitude-longitude and latitudealtitude. In each grid, the Es ORs are obtained by dividing the number of the Es-related S4max values by the number of the qualified S4max values. In order to reduce the random errors in deriving the Es ORs, it is set that for each grid, only when the number of available qualified S4max values is greater than 5, the Es ORs will be calculated using the equation shown in Fig. 3, otherwise, the grid will be left blank. What needs to be mentioned is that in previous studies, although the distribution patterns of Es ORs derived from GPS RO data are generally in accordance with each other, the magnitudes of Es ORs vary significantly. The global Es ORs derived by Arras et al. (2008) based on the 50-Hz signal-noise-ratio (SNR) data vary between 0 and 50%, while those derived by Chu et al. (2014) using the 1 Hz and phase data vary between 0 and 30%. In the work of Tsai et al. (2018), which uses the 50 Hz SNR and phase data, the Es ORs over some regions are greater than 80%. Differences in the RO data types, in the grid size and in the detailed methods used to identify Es layers would lead to the differences in the magnitudes of the derived Es ORs, but they would not influence the general features of the climatology (Niu 2021). Figure 4 presents the monthly variations of the global latitude-longitude distribution of the Es ORs during 2007 to 2015 binned with 2° × 2° grids. The magnetic inclination contours, which are calculated using the IGRF12 model, are also shown in each panel. It can be seen from Fig. 4a-l that the global distributions of the Es ORs show significant monthly variations and the maximum (minimum) occur in the middle latitude regions in local summer (winter). In the northern hemisphere (NH) summer, for the middle latitude regions, the Es ORs are the highest over the Southeast Asia and the lowest in the North America; while in the southern hemisphere (SH) summer, the Es ORs are the lowest over the southern Atlantic anatomy (SAA) zone compared with other regions in the same latitude. For the low latitudes, seasonal dependency of Es ORs is weaker than at middle latitudes, especially, the Es ORs are generally low around the magnetic equator in all months. For the high latitudes, the Es ORs are quite low compared with middle latitudes, and due to the uneven data distribution and our threshold limit, there are a few blank spaces in the north and south poles.

Fig. 3 Flowchart of Es ORs calculation
The red curve in each subfigure of Fig. 4 represents the contour line with a magnetic inclination of 0. It can be seen that in each month, along the red curve, there is a reduction in Es ORs inside a narrow zone of about 3-5 degrees in latitude centered at the equator. Besides, there are significant reductions of Es ORs at high latitudes where inclination angles are larger than 70-80 degrees. These phenomena indicate that due to the unique magnetic field geometry and the need for the plasma to remain neutral, too large or too small a magnetic dip will inhibit the formation of the Es layers (Haldoupis et al. 2011). Figure 5 depicts the magnetic latitude (mag-latitude)time variations of the Es ORs binned with a grid of 1° × 5 days during the period from January 2007 to December 2015. Based on the grids with high spatial and temporal resolutions, the expected summer maximum alternating between the northern hemisphere and the southern one is clearly visible. It can be seen that the Es layers mainly occur at middle and low mag-latitudes between 15°N and 55°N of the NH and between 20°S and 60°S of the SH in local summer seasons, which is consistent with previous studies (Arras et al. 2008;Tsai et al. 2018;Yu et al. 2019). It can be seen from Fig. 5 that Es ORs are much higher in the two extended middle mag-latitude zones than that in the magnetic equatorial zone, and are lowest in the two auroral zones. Besides, the hemispherical asymmetry of the Es ORs is apparently shown in Fig. 5, i.e., the Es ORs at NH middle mag-latitudes during NH summers are greater than those at SH middle mag-latitudes during SH summers. This hemispheric asymmetry is also found by Arras (2010), Tsai et al. (2018) and Li et al. (2020), which is mainly due to the difference in the strengths of the horizontal magnetic field in the two hemispheres. More specifically, the distribution characteristics of the horizontal magnetic field leads to the low value of Es ORs in the Southern Atlantic region and surrounding areas. Since Es ORs are integrated over all longitudes in each 1° mag-latitude range to reach the values shown in Fig. 5, the low Es ORs over the extended SAA diminish the integrated Es ORs over the corresponding mag-latitude zones significantly.
In order to show more clearly the intensity of Es ORs over time, we have used a 30-day moving average method to calculate and to plot the time-varying curves of Es ORs at different mag-latitudes, as shown in Fig. 6a, b. A line plot of the average Es ORs with mag-latitude are also shown in Fig. 6c for better understanding the statistical relationship between mag-latitude and Es ORs. In Fig. 6a, b, the blue and the green curves represent the Es ORs in middle mag-latitude regions, which have the largest values compared to those in other latitudinal regions and have distinct seasonal variations. We found that in the NH (Fig. 6a), the maximum of Es ORs reaches about 80%, while in the SH (Fig. 6b), the maximum of Es ORs reaches about 70%. In both of the two subfigures, the red curves represent the temporal variations of the Es ORs in the magnetic equatorial region. It can be seen that the Es ORs in this region are significantly lower than those in the middle and low mag-latitude regions, and the seasonal variations are not significant. The gray curves represent the temporal variations of Es ORs in the high mag-latitude areas. It can be found that compared with the magnetic equatorial area, the Es ORs in high mag-latitude areas are generally of lower magnitudes. Figure 6c counts the average value of Es ORs for each degree of mag-latitude in the range of 85°S-85°N. The results show that in the geomagnetic equatorial region, the average value of Es ORs is reduced compared to the mid-latitude and low-latitude regions. The smallest average of Es ORs occurs in the polar regions, but there is an increase in the area around 75°N, which is consistent with the distribution of Es ORs shown in Figs. 4, 5. In addition, compared with Fig. 5, Fig. 6c seems to have a latitude drift, which is actually a deception. This is caused by averaging the summer maximum and winter minimum at the same mag-latitude in Fig. 5. For the 85-90° mag-latitude areas, the results may lack reliability due to the small amount of observational data, so they are not discussed here. Figure 6a, b also shows that the Es ORs in middle maglatitude regions of the NH are generally significantly higher than those in the corresponding mag-latitude zone of the SH. The Es ORs also have inter-annual variabilities. The maximum values of the Es ORs shown in the blue and the green curves are smaller during 2013-2015 compared with those during 2009-2012. The physics behind Es layers relies on the "wind shear" theory which is now widely accepted as the mechanism responsible for Es layer formation (Whitehead 1989). According to this theory, vertical shears in the horizontal wind can drive the long-lived metal ions in the lower thermosphere to move vertically and converge into dense plasma layers by the combined action of ion-neutral collisional coupling and magnetic Lorentz forcing. It is expected that the factors affecting the formation of Es layers mainly include neutral wind, magnetic fields, electric fields and metal ion concentration. Therefore, the possible differences in factors such as wind field, magnetic field, and meteor flux should be considered as the cause of the difference in Es ORs between NH and SH regions. Since the average annual variation of the magnetic field is insignificant, it can be excluded as the cause of the inter-annual variation of Es ORs. Considering that the atmospheric wave dynamics and the concentration of metal ions play key roles in the formation process of Es layers by influencing the vertical wind shears and by providing substance needed for convergence (Mathews 1998;Haldoupis and Pancheva 2006;Haldoupis et al. 2007), the anomalies of the wind field in the atmosphere (Chang et al. 2018) and the concentration of metal ions should be the factors responsible for the inter-annual variation of Es ORs, and further investigations on these topics are needed in the future.
To further illustrate the influence of the magnetic field strength on the formation of Es, Fig. 7 is provided to show the global distribution of the intensity of horizontal magnetic field (B h ) in the left panel and the magnetic field inclination in the right panel, computed from IGRF12 at 105 km altitude. It is visible in Fig. 7a that large B h values are located in low and middle latitudes with a strong maximum over southeast Asia, and the horizontal intensity decreases quickly poleward. A comparison of Fig. 7(a) and Fig. 4 implies that the B h is one of the reasons behind the distribution of Es ORs. As pointed out above, Es ORs increase at middle latitudes during summer seasons with two noticeable exceptions of minima zones, one over the SAA zone in the SH and the other over northern America in the NH. In these two regions, the B h is strongly reduced, and the magnetic inclination also shows atypically high values in SAA area (cf. Fig. 7b). The global map of the B h indicates a decisive role of horizontal magnetic field in Es formation through the Lorentz force which drives the ion convergence in the zonal wind shear mechanism. In addition, from Fig. 7a we can see that the B h in the NH is significantly greater than that in the SH. Therefore, it can be inferred that the difference in B h is one of the fundamental reasons for the difference in Es ORs between the northern and southern hemispheres. In Fig. 7b, it can be seen that the inclination of the magnetic field in the equator area is close to 0, so that the electrons are confined to magnetic field lines that are exactly horizontal along the equator. Therefore, they are not able to follow converging ions moving with neutral winds. As a consequence, a large electric field would build up prohibiting the Es formation process (Kelley et al. 1989;Arras 2010), and then the low Es ORs along the geomagnetic equator are explained. Besides, in polar latitudes, the magnetic dip angles are large, thus the B h which is involved in the wind shear process is small (cf. Fig. 7a). As a result of these geometrical constraints, wind shear fails in polar latitudes, where only a small amount of Es layers can be observed (Haldoupis 2011). Figure 8 presents the monthly zonally averaged values of the altitude-magnetic latitude distributions of the Es ORs binned with the grids of 2 km altitude ×2° latitude during the nine years from 2007 to 2015. As shown, the distribution of Es ORs in each month is apparently divided along the two sides of the geomagnetic equator, and most Es layers occur in the middle mag-latitude regions, with the central height range of 100 km to 120 km in the summer hemisphere and 100 km to 115 km in the winter hemisphere. Within these height ranges, the maximum Es ORs of about 80% and 70% appeared during June in the NH and during December in the SH, respectively. The patterns of seasonal variations in the NH and SH are similar. Take the NH as an example, it can be seen from Fig. 8a-f that from spring to summer, the Es ORs in the middle and low mag-latitude regions gradually increase, and the spatial coverage of the Es layers gradually expands, whereas the season variation pattern of the Es ORs from autumn to winter (Fig. 8g-l) is opposite. Figure 9 presents the altitude-local time (LT) distributions of the Es ORs in different mag-latitudes of the two hemispheres for the four different seasons. Here, the seasons are categorized as MAM (Mar., Apr. and May), JJA (Jun., Jul. and Aug.), SON (Sept., Oct. and Nov.) and DJF (Dec., Jan. and Feb.). Taking the NH as an example, from MAM to JJA (Fig. 9a1-b6), the mag-latitude zone of the Es ORs maximum values move from 15 to 30°N (MAM) to 30-45°N (JJA), however from SON to DJF (Fig. 9c1-d6), the minimum values of Es ORs extend from the higher mag-latitudes to the lower ones. The diurnal variations of the Es ORs in different seasons over different mag-latitudes are clearly presented in Fig. 9, which are mainly characterized by two maximums (hereafter called semidiurnal behavior) in the morning and afternoon to evening in the mag-latitudes range of about 15-60° in the summer hemisphere during JJA and DJF, and one maximum (hereafter called diurnal behavior) in the afternoon, which is dominated throughout the year in the low mag-latitudes of 0-15° and in local winter months of other mag-latitudes. More diurnal behaviors are also discernable in the equinox seasons. Figure 9 also shows that the Es layers appear most often in the daytime, in the areas with an altitude range of 100-110 km. In the middle and low mag-latitudes (15-45°) of the summer hemisphere, the areas where Es layers most appear extend to 100-120 km, and the Es layers begin to descend to lower altitudes from 120 km at about 18 LT, and eventually tends to stay at around 100 km, which is basically consistent with recent studies (Shinagawa et al. 2021;Andoh et al. 2021). Andoh et al. (2021 found that the temporal evolutions of 3-D metal ion layer (MIL) structures around Japan can be classified into four phases, and in the last phases, MILs mainly stagnate at 100 km or continue to descend to lower than 100 km, depending on the strength of the downward winds. In most subfigures of Fig. 9, Es ORs basically stagnated around 100 km, which is consistent with the conclusion of Andoh et al. (2021). In addition, Yu et al. (2021b) reported a study of ionospheric irregularities using scintillation data from COSMIC satellites and identified a large-scale horizontal transport of long-lived metallic ions in Es layers. They proposed a "generalized wind shear theory" with threedimensional motions of ion drift, and pointed out that the lower thermospheric meridional circulation influences the meridional transport and seasonal variations of metallic ions in Es layers. In Figs. 8, 9, the latitude zone where the maximum of Es ORs is located changes with months. This is a typical phenomenon of meridional transport of Es layers, and can be explained by the winter-to-summer interhemispheric transport of ions caused by the lower thermospheric meridional circulation. Figure 10 presents the magnetic latitude-LT distributions of the Es ORs at different altitude ranges in the four different seasons. As shown, Es is mainly formed at 100-110 km, and the peak altitude of Es ORs extends to 100-120 km in the middle mag-latitude areas of the summer hemisphere. By carefully comparing the subfigures of Fig. 10, it can be seen that 100 km is the downward boundary of maximum Es ORs. Depending on the strength of the downwind, the Es ORs observed in the altitudes lower than 100 km are significantly reduced, which is consistent with the previous analysis. The characteristics of the diurnal variation of the Es ORs in different seasons at different altitudes are also apparently shown in this figure. Taking the NH summer as an example, as shown in Fig. 10b1, the diurnal behavior is dominant in the low altitude range of 90-95 km, while the semidiurnal behavior can be identified in the altitude ranges of 95-115 km. Figure 11 shows the local time-month variation of the Es ORs at low and middle latitudes during January 2007 to December 2015. We use the sunrise and sunset times at sea level as a simple reference for the beginning of day and night. As shown, in the middle and low latitudes, the Es ORs increase after sunrise and decrease after sunset. In another word, Es mainly occurs during the daytime. Moreover, Es ORs in mid-latitudes have obvious semidiurnal behavior during daytime. The first peak appears before noon (6-10 LT), and the second peak appears around sunset (16-22 LT). The equator and its surrounding areas have a single-peak characteristic, with the peak appearing before sunset (16-18 LT). During night, in summer of mid-latitudes (cf. Fig. 11a, d), the average value of Es ORs is still greater than 50%. In winter of midlatitudes, the Es layers are rarely observed at night. In the equator and its surrounding area (cf. Fig. 11h), Es ORs have the lowest values from midnight to sunrise, being about 25%. Figure 12 shows the local time-month variations of the Es ORs at polar regions. As shown, the Es ORs during the polar day are generally greater than in polar night. Taking Fig. 12a, b as an example, in polar day, the maximum value of Es ORs is about 60%, and most of the values are higher than 30%; in polar night, the maximum value of Es ORs is about 50%, and most of the values are below 30%. Similarly, in other regions shown in Fig. 12, when these regions are in polar day, the range of Es ORs is 30-50%. However, when these regions are in polar night, the values of Es ORs In each panel (a-i), the red and the black curve represents the sunrise and sunset time in each month for the specific latitude zone are generally lower than 30%. An interesting finding is that in the 75°N-80°N latitude areas (Fig. 12a, b), the Es ORs increase between 18 and 24 LT, which explains the abnormal increase of Es ORs in the 75°N area shown in Fig. 6c, but its physical mechanism remains to be studied.

Wind shear simulation
The current mainstream view on the formation of the mid-latitude Es layers relies on the wind shear mechanism. For the formation of Es layers, charged particles including metal ions are compressed into a thin layer by vertical wind shears. In this process, metal ions are mainly driven by horizontal winds, ion-neutral collision and the Lorentz force of magnetic field, and then electrons follow them moving along the magnetic field lines to maintain charge neutrality. One of the unsolved issues in wind shear simulation is that it does not well explain the overall distribution and seasonal variation of Es layers. In this section, in order to further explain the spatial distribution and temporal variation of the Es layers, we studied the monthly convergence of vertical ion velocity of 100-115 km. This height range is the intersection of the center heights of Es in summer and winter seasons, as shown in Fig. 8, and was found by Yu et al. (2019) to be feasible for wind shear simulation. The vertical ion drift velocity caused by vertical shear of the horizontal wind of the neutral atmosphere can be written as follows (Mathews 1998): (1) w j = V j cos I j sin I j + r j U j cos I j 1 + r 2 j , where U j and V j represent the eastward neutral velocity and the southward neutral velocity of horizontal wind components, respectively, which were calculated by HWM14 wind field model in the present study. I j represents the magnetic dip angle, which can be calculated by the IGRF12 model. r j represents the ratio of the ion-neutral collision frequency to the ion gyrofrequency. In order to simplify the calculation, we derive the r j in j region based on the altitude profile of r 0 provided by Bishop and Earle (2003), as shown in formula (2): where B mean is the average value of the total magnetic field strength, and B j is the total magnetic field strength in the j area. Note that r j is a rough estimate, and the effects of ionospheric electric field, the vertical motions of gravity waves and the meteor influx rate are neglected (MacDougall et al. 2000a, b;Yeh et al. 2014;Yu et al. 2015). The vertical ion convergence (VIC) is defined as follows: Arras (2010) and Yu et al. (2019) indicated that sporadic E layers appear preferably in regions of negative vertical zonal wind shear ( ∂w i ∂z < 0 ). Therefore, in combination with formula (3), only positive values of the VIC contribute to Es layer formation (Shinagawa et al. 2017;Qiu et al. 2019). In our analysis, VIC during . Fig. 12 The same as Fig. 11, but for the polar regions. In each panel (a-f), the areas enclosed by the red lines represent the polar day, and the areas enclosed by the black lines represent the polar night the period from January 2007 to December 2015 was calculated by Eq.
(3). The calculated VIC is binned and averaged in a 2° latitude ×2° longitude grid at altitude range from 100 to 115 km. Figure 13a-l presents the global distributions of the simulated monthly variation of VIC in units of 10 −3 s −1 , which are of good correlations with the geographical distributions of the Es ORs derived from COSMIC RO S4max data shown in Fig. 4. VIC is mainly distributed in the middle and low latitudes, while in the North America and the SAA region, VIC is weakened. VIC generally increases in local summer and weakens in local winter. In spring and autumn months, VIC is distributed on both sides of the equator and extends to high latitudes. The monthly variations of the distributions of the simulated VIC values support the wind shear theory of Es formation and demonstrate that the seasonal dependence and the overall distribution of Es layers are partly attributed to the convergence of vertical ion velocity driven by neutral wind.
However, we noticed that some disagreements exist between the latitude distributions and the seasonal variations of VIC and those of Es ORs. VIC generally extends to the equator and high latitudes, occupying a wider range of latitudes than Es ORs. In addition, there are also differences between the monthly variations in VIC and the Es ORs in the equinox months. Yeh et al. (2014) used the HWM07 model to simulate the wind shear theory and explained the formation of the Es layers in combination with the mechanism of meteor ionization. They pointed out that in DJF months and JJA months, the distribution of Es layer activity is correlated with the distribution of the convergence mechanism, which is associated with wind shear theory. But in equinox months, they found that the correlations are not obvious, and cannot be explained using wind shear theory along but can be explained using the meteor ionization mechanism. Figure 14 shows the annual variation of VIC with magnetic latitudes. The results show that the simulated VIC has similar variation characteristic to Es ORs (cf. Fig. 5), VIC and Es ORs all increase in summer at middle and low mag-latitudes. However, the distribution of VIC at the geomagnetic equator is different from the distribution of Es ORs. This is because the electrons in this area are confined to magnetic field lines that are exactly horizontal along the equator, which gives them resistance when they follow the converging ions moving with Fig. 13 Simulation results of monthly variations of the global distribution of vertical ion convergence derived from the HWM14 and IGRF12 models at altitudes ranging between 100 and 115 km, during the period from January 2007 to December 2015 neutral winds (Kelley et al. 1989;Arras 2010). In addition, Fig. 14 also shows that the VIC simulated based on the HWM14 model cannot provide accurate annual variation information. The main reasons are that the model does not consider the annual variation of metal ion concentration (Yeh et al. 2014;Chu et al. 2014) and the HWM14 model fails to provide sufficiently detailed information about the annual variation of the wind field (Drob et al. 2015).
We already know that the semidiurnal behavior of Es ORs mainly occur in summer in mid-latitudes (cf. Figs. 9 and 10). Figure 15a shows an example of the diurnal variation of VIC in summer in middle mag-latitudes (30-45°N) in the NH. As shown, VIC has obvious semidiurnal behavior, with the first peak occurring before noon, and the second peak occurring in the afternoon to evening. Figure 15b shows another example of the diurnal variation of VIC, and the results also show that VIC has a semidiurnal behavior in summer in middle mag-latitudes. The diurnal characteristics of VIC are basically consistent with the distribution of Es ORs (cf. Figs. 9b3, 10b5 and 11a, b). Shinagawa et al. (2017) studied the altitude profiles of the zonal wind amplitude, and they pointed out that the semidiurnal components of the zonal wind was mainly observed above 110 km, and the zonal wind peaks appear in the morning and afternoon to evening, respectively. Therefore, for the semidiurnal behavior of Es ORs, a reasonable explanation is that the zonal wind strengthens in the morning and afternoon to evening, leading to an increase in VIC, and finally making Es ORs peak during these two time periods. In addition, Fig. 15a also shows that at an altitude of 100 km, the VIC becomes very small, which explains that the Es layers are rarely observed at an altitude of 100 km. It should be noted that since the HWM14 model uses the universal time for calculation, Fig. 15 has converted the universal time to the local time.

Conclusions
In this study, the global distributions and seasonal variations of Es ORs are investigated by using the S4max data from COSMIC RO measurements during 2007 to 2015. The latitudinal and seasonal dependences of Es ORs are identified. The detailed characteristics in the spatial and temporal distributions of Es ORs are presented, including the summer maximum in Es ORs at middle latitudes, the low Es ORs in the North America and the SAA region, the infrequent occurrence of Es layers over the geomagnetic equator and the minimum Es ORs in the polar region. In addition, an abnormal increase in Es ORs was observed in the 75°N area, which has not been reported before. The present study analyzes the relationship between Es ORs and the Earth's magnetic field, and the results show that the Earth's magnetic field, including its horizontal intensity and inclination, is one of the factors responsible for the global distribution of Es layers. In areas with high B h , Es ORs usually increase. However, in the equatorial region, the inclination of the magnetic field is close to 0, which prevents the compression of electrons or ions near this area, resulting in a decrease in Es ORs in this area. In addition, the experiment also shows that the average horizontal magnetic field strength in the NH is greater than that in the SH, which is one of the reasons for the difference in Es ORs between the two hemispheres. The analysis of the diurnal variation of the Es ORs shows that the Es layers mainly occur during the daytime. In particular, we also separately analyzed the Es layers in the polar regions. Experiments show that Es ORs increases during the polar day and decreases during the polar night. The prominent characteristics of the diurnal variations of the Es layers also include the diurnal behavior and semidiurnal behavior of the Es layers. The semidiurnal behavior of the Es layers is common in summer in middle mag-latitude regions. As time goes by, the peak altitude of the Es layers drops from about 120 km to about 100 km. In equinox and winter months, the diurnal behavior of Es ORs occurs more common than the semidiurnal behavior over the globe. To understand the mechanism underlying the spatial and temporal variations of Es layers, the VIC of Fe + are simulated based on wind shear theory. It is found that the distributions of the simulated VIC induced by the neutral wind shear at the altitude range of 100-115 km could partially explain the latitudinal dependence and the summer maximum (or winter minimum) of Es ORs. In addition, the analysis of the diurnal variation of VIC shows that VIC has semidiurnal behavior in summer over mid-latitudes, and its peaks are mainly distributed above 110 km altitude. But the value of VIC is significantly reduced at 100 km altitude. The experiment shows that the diurnal variation characteristics of VIC can explain the semidiurnal behavior and height distribution of Es layers in summer over mid-latitudes. The disagreements between the variations of VIC and Es ORs indicate that other processes, such as the meteor influx rate, ionospheric electric fields and atmospheric tides, should also contribute to the formation and the variations of Es ORs.
The present study is expected to contribute to the understanding of the global distribution, the seasonal variation and the influencing factors of the Es layers, and it provides positive support for the current mainstream views about the formation mechanism of Es layers. However, the variabilities in the Es layers caused by the concentration of metal ions and atmospheric tides need further investigation other than a brief mention. In addition, the present study hasn't provided explanations for the abnormal increase of Es ORs in the area of approximately 75°N, so further investigation is needed for the formation of the Es layers in high latitudes.