Introduction

Traveling ionospheric disturbances (TIDs) are known as the manifestation of atmospheric gravity waves (AGWs) of lower atmospheric origins (Hines 1960) in the ionosphere. Among the TIDs, the disturbances with horizontal wavelength of 100–300 km and phase velocity of 50–300 m/s are classified as medium-scale TIDs (MSTIDs) (e.g. Hunsucker 1982; Otsuka et al. 2013). MSTIDs were previously studied using various techniques in the middle latitudes such as all-sky airglow imagers, ionosondes, incoherent scatter radar, and total electron content (TEC) observations derived from ground-based global navigation satellite system (GNSS) receivers. (e.g., Bowman 1985, 1989; Garcia et al. 2000; Saito et al. 2002; Shiokawa et al. 2003a, 2006; Otsuka et al. 2004, 2011; Kotake et al. 2007; Tsugawa et al. 2007b; Amorim et al. 2011; Lakshmi Narayanan et al. 2014; Medvedev et al. 2017; Takeo et al. 2017).

MSTIDs were previously categorized into daytime and nighttime types according to their wave characteristics. The daytime MSTIDs at mid-latitude region in northern hemisphere have feature of frontal alignment in east–west (E–W) direction propagating in meridional north–south (N–S) direction or frontal alignment of southwest–northeast (SW–NE) direction propagating in northwest–southeast (NW–SE) direction (Hernandez-Pajares et al. 2012; Kotake et al. 2006; Otsuka et al. 2011). On the other hand, the nighttime MSTID has a unique frontal alignment in the northwest–southeast (NW–SE) direction with propagation direction in equator westward direction (Behnke 1979; Saito et al. 1998; Garcia et al. 2000; Shiokawa et al. 2003a; Kotake et al. 2006; Tsugawa et al. 2007a; Kubota et al. 2011; Hernández‐Pajares et al. 2012; Rajesh et al. 2016). It is due to such a special frontal alignment, the phenomenon was originally explained by Perkins instability (c.f. Perkins 1973; Behnke 1979; Garcia et al. 2000).

As the growth rate of Perkins instability is much slower than the observations, previous studies suggested that the instability also requires seeding perturbation to induce the instability development. Kelley and Fukao (1991) recommended that AGWs could be a seeding mechanism to accelerate the Perkins instability. Huang et al. (1994) suggested that AGWs play an important role in formation of the nighttime MSTIDs due to the electrodynamic coupling, which was supported by follow-up observations (Miller et al. 1997; Nicolls and Kelley 2005). Meanwhile, Shiokawa et al. (2003b) calculated the theoretical predicted growth of the Perkins instability. The result indicated that the Perkins instability has low growth rate alone and it takes more than 2 h, which is much longer than the observations. Recently, Chou et al. (2017, 2018), employed dense GNSS observations and simulations on the study of MSTIDs driven by typhoon and proposed that the severe weather-driven concentric GWs (CGWs) could accelerate the Perkins instability by providing additional polarization electric field perturbations driven by the perturbation winds of CGWs.

Additionally, Sporadic E layer (EsL) instability is proposed as an important indirect driver of the nighttime MSTIDs (Cosgrove and Tsunoda 2004), which is mainly produced by the effect of vertical wind shear in E region (Whitehead 1961; Mathews 1998; Carrasco et al. 2007; Haldoupis 2011; Yeh et al. 2014). EsL instability can act as a vital character for accelerating growth of Perkins instability by projecting the polarization electric field perturbations to the F region. Tsunoda (2006) showed the analytic results on the acceleration of the instability and Yokoyama et al. (2009) showed the self-consistent coupling of EsL and F layer Perkins instabilities through numerical simulations. Both studies have indicated the importance of this E and F layer coupling for growth of the MSTIDs.

Not only occurring at mid-latitude, observational studies also depicted characteristics of MSTIDs over the low-latitude regions (c.f. Shiokawa et al. 2006: Lee et al. 2008; MacDougall et al. 2011; Fukushima et al. 2012; Lakshmi Narayanan et al. 2014; Jonah et al. 2016; Paulino et al. 2016; Figueiredo et al. 2018). For instance, the investigations using OI 630.0 nm airglow images taken by all-sky imager around the transition region of low to mid-latitudes (19.3° N dip latitude) reported by Lakshmi Narayanan et al. (2014) suggested that the MSTID could be inhibited by the equatorial ionospheric anomaly. Shiokawa et al. (2006) and Fukushima et al. (2012) suggested that the directionalities of MSTIDs over Kototabang (10.4° S dip latitude) are different from the typical MSTIDs observed at the mid-latitudes and attribute their formation by gravity waves. Paulino et al. (2016) carried out an investigation at the off-equatorial region, over São João do Cariri, using airglow observations over nearly a solar cycle and speculated that the gravity waves are the likely source of the periodic MSTIDs.

As described above, MSTIDs were statistically studied in the ionosphere for quite a long time in mid-latitude region with some additional documentations over low-latitude region possibly driven by various physical mechanisms. Most of the studies, however, were performed mainly by adopting visual inspection of the observation data that requires time consuming efforts and could be influenced by subjective bias while inspecting the airglow images or TEC maps. Takeo et al. (2017) performed three-dimensional fast Fourier transform (3-D FFT) to detect the MSTIDs automatically in airglow imagers over a solar cycle for the first time. Here, our main goal is to develop an autonomous detection algorithm by combining the 3-D FFT and support vector machine (SVM) to distinguish and categorize MSTIDs from other waves or irregularities in the GNSS-TEC maps.

In this study, the autonomously identified MSTIDs at low to mid-latitudes are statistically organized according to their wave parameters, i.e., wavelength, amplitude, propagation direction, seasonal and local time (LT) dependences. The possible mechanisms responsible for these low-latitude MSTIDs are discussed based on these statistical results with aids from the EsL observation obtained from the radio occultation (RO) observations of FORMOSAT-3/COSMIC (Yeh et al. 2014) and the AGWs simulated by the high-resolution Whole Atmosphere Community Climate Model (WACCM) (Liu et al. 2014).

Methodology

We develop an autonomous algorithm to detect MSTIDs by applying both the 3-D FFT and the SVM on the GNSS data derived from ~ 150 receivers in Taiwan. As the equatorial plasma bubbles (EPBs) affected by zonal neutral winds could appear in the reverse C-shape spatial distribution (Shiokawa et al. 2004; Kil et al. 2009) showing NW–SE (NE–SW) waveform in the northern (southern) hemisphere similar to the frontal alignment of MSTIDs, it is important to have the algorithm capable of distinguishing the two different instabilities. The algorithm is applied on the GNSS data around Taiwan (20–30° N, 115–125° E geographic latitude/longitude and 15° N dip latitude) under the northern EIA crest during 2013–2015 to study the characteristics, e.g., propagation direction, LT and seasonal dependence, of MSTIDs and EPBs.

Before carrying out the 3-D FFT and SVM, a 60 min. high-pass filter is applied to the GNSS-TEC data followed by interpolation of the filtered TECs into the fixed grids with spatial resolution of ~ 0.5°. latitude/longitude in both zonal and meridional direction with an accumulation period of 30 min. It should be noted that we substitute zero to the region with no available data owing to the limited coverage of GNSS observation around Taiwan. The ionospheric pierce point (IPP) for converting slant to vertical TEC is set at 300 km altitude. An example of MSTID detected by the autonomous algorithm developed in this study is presented in Fig. 1a using filtered TEC over the region, and an example of EPB detection is presented by utilizing rate of TEC index (ROTI) in Fig. 1b. Figure 1c illustrates the grids for the autonomous detection algorithm.

Fig. 1
figure 1

Example of MSTIDs and EPBs detected by the autonomous algorithm developed in this study and are, respectively, presented in a a 60 min. high-pass filtered GNSS-TEC for highlighting the features of MSTIDs and b the rate of TEC index (ROTI) for highlighting the structure of EPBs. The horizontal ranges and grids of the ionospheric phenomena detection (red lines) and the location of the GNSS stations (blue triangles) used in this study are presented in c with the latitude and longitude range of 20–30° N and 115–125° E and ~ 0.5 degree resolution for both zonal and meridional directions

Full size image

We exploit the 3-D FFT method similar to that developed by Matsuda et al. (2014) to the cumulative three-dimensional filtered TEC data for the zonal and meridional wave numbers and angular frequency. It is noted that positive wave number represents northward or eastward direction of propagation. The phase velocities (\(V_{{\text{p}}}\)) of the waves are given in terms of the parameters obtained from the 3-D FFT as follows:

$$V_{{\text{p}}} = \sqrt {\left( {\frac{\omega }{2\pi }\frac{{W_{{\text{m}}} }}{K}} \right)^{2} + \left( {\frac{\omega }{2\pi }\frac{{W_{{\text{z}}} }}{L}} \right)^{2} } ,$$
(1)

where Wm and Wz are meridional and zonal width of the detection region; K, L and \(\omega\) represent meridional wave numbers, zonal wave numbers and angular frequency, respectively. The propagation direction (\({\uptheta }\)) is given as

$$\theta = \emptyset + \tan^{ - 1} \frac{{\left( {\frac{{W_{{\text{m}}} }}{L}} \right)^{2} }}{{\left( {\frac{{W_{{\text{z}}} }}{K}} \right)^{2} }},$$
(2)

where \(\emptyset\) is the azimuth angle, which depends on the value of both zonal and meridional wave numbers.

The 3-D FFT will provide a weaker artificial enhancement of power spectral density (PSD) in the opposite direction than that of the true direction (Takeo et al. 2017). Therefore, the maximum PSD is taken to determine the direction for each computation to rule out the wrong directions in this study.

From the computed results together with the TEC observations, we obtain ten kinds of parameters including horizontal wavelength, meridional and zonal wave numbers, phase velocity, propagation direction, angular frequency, PSD of the 3-D FFT, minimum TEC, day of year (DOY) and universal time (UT). These parameters are exploited for training the SVM model to distinguish MSTID, EPB and quiet time background ionosphere. The detailed procedures for creating the SVM model is illustrated in Fig. 2. In the SVM model, we create a ten-dimensional space where the coordinates stem from the ten kinds of parameters mentioned above. The data are represented as points and are mapped into the ten-dimensional (10-D) space according to the corresponding values of each parameter. In the 10-D space, the SVM model separates different classes by computing the best boundaries called hyperplanes. New data are then mapped into that same space and are predicted into various classes based on the side of the hyperplane on which they might fall. Presumably, the hyperplane between anisotropy classes should be as wide as it could be to maintain the quality of the SVM model. Therefore, a soft-margin method, defining the width of a hyperplane in non-linear classification, is applied on supporting the hyperplane. We define a hyperplane with a soft margin in the following form (Yates 1937; Burges 1998; Chen et al. 2004; Chang and Lin 2011),

$$\left| {N_{{\text{h}}} \cdot x - b} \right| \ge 1 \pm \varepsilon_{i} ,$$
(3)

where \(N_{{\text{h}}}\) represents the normal vector to the hyperplane; x refers to the normalized 10-D parameters; b is a normal constant; the normal constant (on the right-hand side of the equation) “1” or “− 1”, considering the absolute value, indicates two different classes; \(\varepsilon_{i}\) represents slack variable which determine the width of the soft margin. It should be apparent that the greater the value of slack variable is, the larger the breadth of the soft margin would be. To maintain the quality of classification, the slack variable should as small as it could be where we set the value as 0.4 in our case.

Fig. 2
figure 2

Flowchart of the procedure of creating an SVM model for classifying different ionospheric conditions

Full size image

In this study, a supervised learning method is used for classifying the different classes. Supervised learning indicates that we have to label the classes for the raw data to create input–output pairs before training the SVM model. In our case, we label the different classes using discriminants and a little visual inspection for the later validation of the trained results. During the training processes, a solar maximum year, 2000, and a solar minimum year, 2009, are used for training and testing the quality of learning model since MSTIDs and EPBs have an opposite solar activity dependence. Due to the 30 min time resolution in the 3D-FFT computation, the total number of data for these 2 years are around 35,000. The dataset with known labeled classes is randomly separated into two groups including 80% of training data (~ 28,000 piece) and the other 20% of testing data (~ 7000 piece) that were not included in the training process. The training data are used for training the prediction model by computing a hyperplane for categorizing anisotropy classes. On the other hand, the testing data applied on testing the trained model are utilized for calculating the accuracy rate of the SVM model prediction. Here, the accuracy is the ratio of number of correct predictions to the total number of input samples (i.e., the entire 20% testing data). During the testing processes, the original dataset of 35,000 total data number is repeatedly split into the 80% training data and 20% testing data for five times to consecutively test the model prediction and to calculate the accuracy. The final testing results after this process indicate that our model has over 96% accuracy rate for distinguishing MSTID, EPB and ionosphere without irregularity.

Results of autonomous algorithm

In the analysis, we define the four seasons as follows: (1) Spring: from 6 February to 6 May for 3 months, centering on the spring equinox; (2) Summer: from 6 May to 6 August for 3 months, centering on the summer solstice; (3) Fall: from 6 August to 6 November for 3 months, centering on the fall equinox; and (4) Winter: from 6 November to 6 February for 3 months, centering on the winter solstice.

Since the EPB is a well-known phenomenon in the low-latitude ionosphere which has been investigated intensively, it is conceivable to validate the characteristics of EPBs using our algorithm, such as occurrence rate, propagation direction and seasonal dependence. Here, the occurrence rate is defined as the ratio between the number of events in which EPBs are observed (M) and the amount of data in which GNSS-TEC exist (N). Since the GNSS-TEC data are available every 30 s during 2013–2015, the value of N is identical in each computation of occurrence rate.

The observed characteristics of EPBs during 2013–2015 can be summarized as follows (figure not shown). First, the EPBs are usually observed at the period of 2000–0200 LT with an eastward propagation direction. Second, the occurrence rate of EPBs is higher in the solar maximum year of 2014, showing positive dependence on the solar activity. Third, the cumulative occurrence rate for these years indicates that the seasonal occurrence rate of EPBs is descending in the order from Spring, Fall, Summer, to Winter. The characteristics of EPBs in their diurnal and seasonal variations are consistent with those reported in the previous studies (e.g. Farley et al. 1986; Sahai et al. 2000; Huang et al. 2002; Pimenta et al. 2003; Fejer et al. 2005; Gentile et al. 2006; Yao and Makela 2007; Paulino et al. 2011; Sharma et al. 2014; Yokoyama et al. 2014; Sun et al. 2016) suggesting the reliability of our algorithm.

Results of MSTIDs

In this section, we display the annual variations (2013–2015) of the characteristics of MSTIDs by excluding the possible contamination coming from the EPBs. The wavelength and horizontal phase velocity characteristics of MSTIDs over Taiwan are, respectively, among 100–400 km and 100–250 m/s, which are consistent with the previous studies (e.g., Hunsucker 1982; Otsuka et al. 2013). The occurrence rate of MSTIDs is calculated using the same approach applied for calculation of the occurrence of EPBs described in “Results of autonomous algorithm” section. The results illustrate that the occurrence and propagation characteristics of MSTIDs have clear seasonal and LT dependences. The occurrence rate of southward MSTIDs (Fig. 3a–c) indicates that they are more discernable during 0800–2100 LT (red line on the x-axis) in Spring and Winter. During Summer, southward MSTIDs mainly occur around 2100–0300 LT. On the contrary, the occurrence rate of northward MSTIDs (Fig. 3d–f) demonstrates that they are more frequently observed around 1200–2100 LT (red line on the x-axis) from Spring to Fall with a secondary peak occurring around 0000–0300 LT (magenta line on the x-axis) between DOY 100–140. In Fig. 3g–i, a transition boundary between the northward and southwestward MSTIDs takes place 3 h after sunset at around 2100 LT in Summer. It is noteworthy that the LT dependence of the southwestward MSTIDs takes place much later after sunset than those reported in the previous literature (Behnke 1979; Saito et al. 1998; Garcia et al. 2000; Shiokawa et al. 2003a; Kotake et al. 2006; Tsugawa et al. 2007a; Kubota et al. 2011; Hernández-Pajares et al. 2012; Otsuka et al. 2013), suggesting an alternative driving mechanism should be proposed. Also, the propagation direction of the northward MSTIDs (Fig. 4) during nighttime illustrates that the MSTIDs mainly propagate northwestward (wavefront of NE–SW direction) and exactly northward (wavefront of E–W direction) which suggests a rational generation mechanism should be raised as well.

Fig. 3
figure 3

Local time and seasonal variations of the MSTID obtained from 3D-FFT and SVM analysis over Taiwan during 2013–2015. It is divided into 10 days and half an hour per grid for calculating the occurrence rate of southward MSTIDs in a 2013, b 2014, c 2015 year and northward MSTIDs in d 2013, e 2014, f 2015 year. The propagation direction of MSTIDs is compared to the dusk terminator (black line) in g 2013, h 2014, i 2015 year

Full size image
Fig. 4
figure 4

The occurrence rate of a northwestward, b northward and c northeastward MSTIDs over Taiwan in nighttime during 2013–2015. The black lines indicate the dusk terminator

Full size image

Discussion

Here we discuss the characteristics of MSTIDs in the low latitude that are unique from those reported in the middle latitude. The three main characteristics at low latitude are briefly described here. First, there are southwestward MSTIDs like those commonly observed in the middle latitude over Japan occurred 3 h after sunset, but the appearance time is a few hours later than those over Japan and the wavelength and amplitude are greater than other MSTIDs over Taiwan. Second, some of the wavefront (or frontal alignment) and occurrence time of MSTIDs are different from those observed in the middle latitude and there are fewer EsL occurring in the low latitude, suggesting necessity of alternative explanations. Third, there are much more occurrences of northward-propagating MSTIDs before 2200 LT and are likely connected to the northward propagation of AGWs. Some insights could be provided by the observation of 557 nm airglow images and the high-resolution WACCM simulation of AGWs. The three main characteristics are discussed as follows.

Southwestward MSTIDs

On the southwestward MSTIDs in the low latitudes appearing ~ 3 h after sunset, AGWs seeding could be a plausible driver as it can accelerate the Perkins instability in low latitude. This effect is particularly prominent when there are concentric AGWs driven by the server meteorological event. Chou et al. (2017) proposed that polarization electric fields driven by the typhoon triggered concentric AGWs could accelerate Perkins instability and confirmed the effect in the follow-on simulations (Chou et al. 2018). With the similar mechanism of AGWs-Perkins coupling, MSTIDs could be generated locally over Taiwan. However, the propagation direction of the MSTIDs should be coherent with the direction of AGWs. In the later section, we will show that northward AGWs are more frequently appeared than the southward ones in Summer and thus an alternative explanation is required.

The alternative explanation is that these southwestward MSTIDs are likely coming mainly from the higher latitude to the lower latitude region over Taiwan instead of being generated locally. This explanation is intuitive as the southward MSTIDs appear later than those observed in the middle latitude. It is also supported by the comparison of the amplitudes and wavelengths of the MSTIDs in the two latitude regions. Comparison of the filtered TEC over Japan and Taiwan in Fig. 5a, b reveals the relationship between the MSTIDs at low and mid-latitudes. To identify the characteristics of the MSTIDs over Japan and Taiwan, we organize filtered TEC as a function of distances and time (Fig. 5c). The filtered TEC along the magenta dashed line in Fig. 5b is chosen for organizing Fig. 5c for studying the characteristics of MSTIDs in this case. The daytime MSTIDs (Fig. 5a) propagate toward southeast and northwest over Japan and Taiwan, respectively. On the contrary, the nighttime MSTIDs (Fig. 5b) demonstrate that they are propagating southwestward over Japan and across a long distance to reach over Taiwan region. In Fig. 5c, it apparently shows that the MSTIDs over Japan have greater wave amplitudes (~ 0.2 TECu, where 1 TECu = 1016  el/m2) and longer wavelengths (~ 250 km) than those over Taiwan where the wave amplitude is around 0.1 TECu and the wavelength is about 100 km. The red rectangle in Fig. 5c, therefore, illustrates that the MSTIDs around 1400–2100 UT (2200–0500 LT) over 22–28° N latitude are likely coming from Japan as they have a greater TEC amplitude fluctuation and wavelengths than those in the same latitudes at earlier local times.

Fig. 5
figure 5

Daytime and nighttime MSTID propagation over Japan and Taiwan. a An example of typical daytime MSTIDs over mid- and low-latitude region, and b example showing nighttime MSTIDs propagating from Japan to Taiwan. c The filtered TEC perturbations, along the purple dashed line in b, organized as a time–distance map

Full size image

To further identify the relationship between the MSTIDs over Japan and Taiwan, Fig. 6 compares the propagation directions with the normalized power spectral density (PSD) of MSTIDs during DOY 100–300 of 2013–2015. The PSD derived from 3-D FFT indicates the power of the waves, which corresponds to the square of the amplitude fluctuation (Danielson and Lanczos 1942; Cooley and Tukey 1965). The normalized PSD with the much greater values, indicating the greater wave amplitudes, mainly occur three hours after sunset (Fig. 6d–f). They exactly match the appearing time of southwestward MSTIDs (Fig. 6a–c). The much greater wave amplitude coincided with the southward MSTIDs suggests that they are very likely originated from Japan and propagating southward instead of being generated in situ over Taiwan.

Fig. 6
figure 6

A comparison between propagation direction and normalized PSD of MSTIDs with a dusk terminator (black line) over Taiwan. The southwestward MSTIDs (deep blue color) in the years a 2013, b 2014, c 2015 is in good agreement with the PSD with greater value (red color) in d 2013, e 2014, f 2015, suggesting that the southwestward MSTIDs have greater amplitudes than those propagating in other directions

Full size image

We suppose that if the southwestward MSTIDs are mainly coming from Japan, the occurrence rate of southwestward MSTIDs over Japan and Taiwan shall be similar. Otsuka et al. (2011) revealed that the occurrence rate is ~ 30% over Japan in a solar maximum year in nighttime in Summer. On the other hand, our study shows a higher occurrence rate of ~ 40% in the year of higher solar activity. The ~ 10% higher occurrence rate may suggest that some of the southward MSTIDs are generated in situ over Taiwan region. It is worthwhile to note that the solar activity condition in our study (f10.7 index is 145.9 s.f.u.) is weaker than the observation period of Otsuka et al. (2011) (f10.7 index is 179.4 s.f.u.), which suggests uncertainties due to the different solar activity level. More analyses comparing the statistics of MSTIDs over Japan and Taiwan during the same period are required in a future study.

Perkins/Es layer instability and MSTIDs around dusk

For the nighttime MSTIDs, Perkins instability (Perkins 1973) was considered as the most plausible generation mechanism to explain the special wavefront alignment (NW–SE) of MSTIDs (Behnke 1979; Garcia et al. 2000). According to Shiokawa et al. (2003a, b) and Tsunoda (2006), the linear growth rate (\(\gamma_{p}\)) of the Perkins instability can be given as

$$\gamma_{{\text{p}}} = \frac{g}{{H_{{\text{n}}} \left\langle {v_{{{\text{in}}}} } \right\rangle }}\frac{{\left. {\sin^{2} \left( I \right)\sin \left( {\theta - \alpha } \right)\sin (\alpha } \right)}}{{\left. {\cos (\theta } \right)}} = \frac{{E_{0} \sin \left( {\theta - \alpha } \right)\sin (\alpha )\left. {\cos (I} \right)}}{{BH_{n} }}\frac{{\left\langle {v_{{{\text{in}}}} } \right\rangle_{{{\text{ot}}}} }}{{\left\langle {v_{{{\text{in}}}} } \right\rangle }},$$
(4)

where g, \(H_{{\text{n}}}\) and \(v_{{{\text{in}}}}\) are the gravitational acceleration, scale height of neutral atmosphere and ion-neutral collision frequency; \(I\) represents the magnetic inclination; \(\theta\) is the angle of background electric field from east; \(\alpha\) means the angle between the direction normal to the frontal structure and east; \(E_{0}\) and B are the background electric field and magnetic field; subscript 0t represents inclusion of a tangential component. Note that \(\alpha\) should be greater than zero but less than or equal to \(\theta\).

In Eq. (4), the ion-neutral collision frequency is inversely proportional to the growth rate. The neutral density is small during solstices and magnetically quiet periods of solar activity, which generates a lower ion-neutral collision frequency. This effect could produce a high occurrence rate of nighttime MSTIDs during solstices and low solar flux years which agree with the previous literature (e.g. Shiokawa et al. 2003a; Kotake et al. 2006).

On the other hand, a coupled electrodynamical system between EsL and F region also plays an important role in accelerating the growth rate for nighttime MSTIDs (e.g., Cosgrove and Tsunoda 2004; Tsunoda 2006; Otsuka et al. 2008; Yokoyama et al. 2009). In this study, we estimate the occurrence of EsL using FORMOSAT-3/COSMIC (F3/C) radio occultation (RO) data for comparison with the MSTIDs characteristics in nighttime during 2013–2015. It is worthwhile to note that such comparison could provide some insights on the coupled system, but it is not necessary to provide the causal link of EsL and MSTIDs occurrence since the wind shear effect is required to produce the EsL instability and the polarization electric fields for accelerating Perkins instability. Such wind shear and whether if the EsL instability actually occurs is not covered in our observations of EsL.

Resende et al. (2018) revealed that the S4 scintillation index derived from the RO profiles of F3/C, which provide the disturbance information from the signal to noise ratio (SNR) of the electron density from GPS L1 band, is suitable to estimate the EsL at low latitudes due to their high correlation. In this case, an S4 threshold of 0.3 is set over 80–130 km altitudes (c.f. Whitehead 1989; Arras et al. 2008, 2009; Chu et al. 2014; Yeh et al. 2014; Tsai et al. 2018) for evaluating the occurrence of EsL. It is noteworthy that EsL does not have a conjugate effect but MSTIDs could be presented as conjugate structures in both hemispheres according to previous studies (Saito et al. 1998; Otsuka et al. 2004; Martinis et al. 2019). That is, the EsL occurring at one of the hemispheres could trigger the MSTIDs in both hemispheres, coherently. As a result, two detection regions are necessary, which are (i) Taiwan region, the same as mentioned in “Methodology” section, and (ii) the conjugate region of Taiwan, the identical meridional and zonal breadth to (i) but centering at 120° E geomagnetic longitude; − 4° S geomagnetic latitude, for investigating the EsL characteristics in nighttime during 2013–2015. It is noted that the occurrences of EsL based on F3/C would be affected by the observation distribution of the radio occultation sounding points as the six microsatellites are continuously moving. Under this condition, the uniform data distribution is important to the results of occurrences. As the six microsatellites of F3/C are separated evenly into six orbital plans performing occultation soundings with GPS, the distribution of the soundings is relatively even comparing to other existing satellite missions. The weakness of our result is that there are fewer soundings available during 2013–2015 comparing with the results from Chu et al. (2014), who used data of 2006–2011 when occultation soundings are the richest in number. Comparing our results of occurrences in the low latitude region to those given by Chu et al. (2014), the general characteristics are similar. Thus our results comparing the relative strength of occurrence rate of EsL along local times and seasons in the low latitude are still useful for comparisons with the occurrences of MSTIDs along local times and seasons.

The EsL over Taiwan and the conjugate region during nighttime (Fig. 7a) show that they mainly occur in Summer with an occurrence rate of ~ 5% where the MSTIDs over Taiwan remain more than 50% of occurrence rate for both south and northward (Fig. 7c, d) with a nearly identical propagation direction before and after sunset (Fig. 7b). However, a clear correlation between EsL and MSTID revealed over Japan in Summer (Otsuka et al. 2008), where the occurrence rate of EsL was greater than Taiwan (c.f. Wu et al. 2005; Chu et al. 2014, Zhou et al. 2017). The nearly absence of EsL together with considerable large occurrences of MSTIDs and their directionality over Taiwan indicate that the EsL-Perkins coupling are less responsible for generating the nighttime MSTIDs over Taiwan in situ. Additionally, the propagation direction of the northward MSTIDs (Fig. 4) illustrates that the MSTIDs mainly propagate northwestward (wavefront of NE–SW direction) and exactly northward (wavefront of E–W direction), which are inconsistent with the theoretical wavefront alignment (NW–SE) of MSTIDs in middle latitude driven by the coupled EsL-Perkins instability (Perkins 1973; Shiokawa et al. 2003a; Cosgrove and Tsunoda 2004). The anisotropy characteristics between EsL and MSTIDs suggest that the coupled EsL-Perkins instability are less likely responsible for generating the MSTIDs around dusk over Taiwan.

Fig. 7
figure 7

Comparison of occurrences of EsL observed by FORMOSAT-3/COSMIC and MSTIDs. a The occurrence rate of EsL detected by FORMOSAT-3/COSMIC RO over Taiwan and the conjugate region of Taiwan in nighttime during 2013–2015. b The propagation direction of MSTIDs over Taiwan in nighttime during 2013–2015. c The occurrence rate of southward MSTIDs over Taiwan in nighttime during 2013–2015. d The occurrence rate of northward MSTIDs over Taiwan in nighttime during 2013–2015. The black lines indicate the dusk terminator

Full size image

Gravity waves and northward MSTIDs

As the characteristics of coupled EsL-Perkins instability do not favor MSTIDs generation and propagation over Taiwan, AGWs are considered as a plausible driver of the northward-propagating MSTIDs in nighttime of Spring and Summer over Taiwan. The secondary peak of occurrence of northward MSTIDs appeared in post-midnight during DOY 100–140 (Fig. 3d–f) should also be amenable to interpretation in terms of AGWs. Figure 8 displays an example of 557 nm airglow images on 29 April 2020 around 0209–0239 LT, obtained from an all sky imager operating at Tainan Astronomical Education Area (120.39° E; 23.18° N). The brightness of the images mainly contributed from about 95–97 km altitude are used to evaluate the horizontal structures of AGWs during this period. Using the 3-D FFT, the series of airglow images demonstrate the salient northeastward propagation of AGWs with a horizontal phase velocity, horizontal wavelength, and period of ~ 114 m/s, ~ 83 km, and ~ 12 min., respectively, over Taiwan. Meanwhile, the filtered TEC map with an IPP at 300 km altitude (Fig. 9) reveals the northeastward MSTIDs with a horizontal phase velocity, horizontal wavelength, and period of ~ 147 m/s, ~ 151 km, and ~ 17 min., respectively. The nearly identical propagation direction and similar characteristics of AGWs and MSTIDs suggesting that AGWs are very likely responsible for the nighttime MSTIDs over Taiwan for this event. However, the time characteristics in the airglow images and filtered TEC maps are different in this case due to limited available observations. Further correlation study between the two observations should be proposed in future work.

Fig. 8
figure 8

The λ557 nm airglow images derived from Tainan Astronomical Education Area on 29 April 2020

Full size image
Fig. 9
figure 9

Two-dimensional filtered GNSS-TEC maps over Taiwan on 29 April 2020

Full size image

To further understand the correlation between AGWs and MSTIDs, we evaluate the resolved AGWs derived from high-resolution WACCM among 30 June to 11 July. As observations of AGWs are very limited in the region of interest, using a sophisticated numerical model that can well reproduce the AGWs could provide some insights for our comparison. Recent studies have shown that high-resolution WACCM simulation could reproduce the AGWs patterns in the mesosphere and lower thermosphere including those semi-circular AGWs triggered by tropical cyclones (Liu et al. 2014; Wu et al. 2018). In the high-resolution WACCM, the horizontal and vertical resolution are ~ 0.25° and ~ 700 m, respectively. As the simulations of high resolution WACCM is expensive in computation, we only compare the observations of MSTIDs with WACCM in Summer when the AGWs effect in triggering MSTIDs is most prominent. In this case, we investigate AGWs horizontal structure at 5.8e−4 hpa (~ 93 km) altitude (c.f. Zhou et al. 2002; Suzuki et al. 2007; Chun and Kim 2008; Heale et al. 2020) in the WACCM. Taking 1400 UT on 9 July for example, the zonal wind (Fig. 10a) and meridional wind (Fig. 10b) patterns display several salient northwestward AGWs structures propagating over north Pacific Ocean and across a long distance ~ 2800 km along the Philippines, before reaching Taiwan. In this case, the waves structures mostly dissipate around 30° N geographic latitude by the stronger in phase background wind. In Fig. 10c, d, we estimate the morphological characteristics of AGWs as a function of time and distances by organizing the zonal and meridional wind along the yellow arrow in Fig. 10a, b. Several long-lasting, ~ 12 h, northwestward AGWs are found during this period with a horizontal wavelength of tens to hundreds of km and a horizontal phase velocity among 20–85 m/s which are consistent with previous literature (Piani et al. 2000; Horinouchi et al. 2002; Sentman et al. 2003; Suzuki et al. 2007; Vadas and Fritts 2009; Vadas and Crowley 2010; Takeo et al. 2017; Tsuchiya et al. 2019).

Fig. 10
figure 10

Examples of a Zonal wind and b Meridional wind derived from high-resolution WACCM at 5.8e−4 hPa at 14:00 UT on 9 July. The yellow arrows indicate the propagation direction of the visible AGWs. c Zonal wind and d Meridional Wind-time-distance map along the yellow arrow in a and b, respectively. Each pair of the red symbols indicates an AGWs event. The numbers above the symbols represent the phase velocity of the AGWs

Full size image

In Fig. 11, the propagation direction characteristic of AGWs (red line with stars), obtained from the 3-D FFT of WACCM results over Taiwan, is compared to the MSTID observations (blue line) in Summer during 2013–2015. The locally generated MSTIDs mainly occurring during 1200–2100 LT in Summer (Fig. 3) are in good agreement with the corresponding WACCM output suggesting that AGWs may play an important role in the generation of MSTIDs during this period. However, there is an apparent disagreement occurring after 2000 LT showing southwestward propagation of the MSTIDs while the AGWs are northward. The disagreement again suggests that the majority of southwestward MSTIDs after 2100 LT mainly come from the southwestward propagation of MSTIDs originated from Japan as described and discussed in “Southwestward MSTIDs” section. It is noted that there is a disagreement between the WACCM results of northward AGWs and the southward turning of MSTIDs around 1000–1200 LT. The northward MSTIDs around 0800–1000 LT turning to southward around 1000–1200 LT are coincided with the much fewer detection of MSTIDs during the period of 0800–1200 LT as clearly seen in Fig. 3g–i. These very few numbers of MSTIDs may be less reliable for comparison. We mark the period of MSTIDs using dashed line to indicate this possible discrepancy. After 1200 LT, the southward MSTIDs are in phase to the southward AGWs until ~ 1500 LT and the results are consistent with those reported by Ding et al. (2011) showing daytime equatorward MSTIDs in the middle latitude.

Fig. 11
figure 11

Comparison between the propagation directions of MSTIDs and the AGWs derived from WACCM. The gray lines are the MSTID propagation direction during 2013–2015, and the blue line represents the averaged value during 2013–2015, while the dashed line indicates the MSTIDs are rarely detected. The red line with stars displays the corresponding WACCM results

Full size image

Tsugawa et al. (2007a) suggested that the MSTID activity over Japan is highest in the nighttime during 2100–0300 JST in Summer (May–August or DOY 140–250). Since the amplitude of the MSTIDs coming from Japan to Taiwan is greater than those locally generated over Taiwan, it is difficult to observe the characteristics of local MSTIDs under the influence of MSTID originated from Japan in the overlapping TEC observations over Taiwan. During DOY 100–140, when the southwestward MSTIDs over Japan were least observed, it could be an opportunity for identifying locally generated MSTIDs over Taiwan. During this time of year, with much fewer MSTIDs coming from Japan, the northward-propagating MSTIDs over Taiwan could be identified and appears as the secondary occurrence peak (Fig. 3). Also, the propagation direction of AGWs is mainly northwestward after 2000 LT (Fig. 11), which is consistent with the secondary peak of the MSTIDs observations (Fig. 7). It suggests that AGWs is likely an important seeding driver for the generation of the northward-propagating MSTIDs over Taiwan in nighttime in Summer. Based on the propagation directions and the seasonal occurrences, we suggest that the AGWs may have an important contribution to the low-latitude MSTIDs in both daytime and nighttime in Summer.

Summary

In summary, we develop an autonomous algorithm that utilizes both the 3-D FFT and SVM to statistically investigate the MSTIDs in the low-latitude equatorial ionization anomaly region over Taiwan during 2013–2015. As EPBs have some characteristics similar to MSTIDs, the occurrence of EPBs is also investigated to distinguish them from the MSTIDs. Several important features such as the variation of the propagation direction and occurrence rates are revealed in this study. The main findings are summarized as follows:

  1. 1.

    The statistical results show that the seasonal, LT variations of TEC perturbations in low latitude over Taiwan region are generally consistent with previous studies of MSTIDs and EPBs, indicating that our autonomous algorithms could successfully distinguish EPBs and MSTIDs in the TEC perturbations.

  2. 2.

    The occurrence rates and propagation directions of MSTIDs show clear seasonal and LT dependences. The southward MSTIDs are observed almost every day around 0800–2100 LT in Spring and Winter; appearing mainly around 2100–0300 LT in Summer; and having the least occurrence in Autumn. On the contrary, northward MSTIDs are observed more frequently around 1200–2100 LT from Spring to Autumn with a secondary peak around 0000–0300 LT during DOY 100–140.

  3. 3.

    The patterns of propagation directions of MSTIDs display a clear boundary of the northward and the southwestward MSTIDs at 2100 LT in Summer. The MSTIDs mainly propagate northward and northwestward around 1200–2100 LT. On the other hand, the MSTIDs observed around 2100–0300 LT mainly propagate southwestward, which might be connected to the MSTIDs originated in the middle latitude region over Japan.

  4. 4.

    By comparing the occurrences, amplitude and wavelength of MSTIDs over Japan and Taiwan, we suggest that the MSTIDs in Taiwan region are influenced by the southwestward-propagating MSTIDs from Japan around 2100–0300 LT during DOY 140–250 due to the overlapping in TEC observations.

  5. 5.

    The wavefront alignment characteristic of northward/northwestward MSTIDs around dusk and within 3 h after sunset is inconsistent with those generated by the coupled EsL-Perkins instability. It suggests that other explanation is required. Comparing the filtered TEC maps with the 557 nm airglow images and simulations of gravity waves by high resolution WACCM illustrates that MSTIDs over Taiwan are likely generated by the AGWs during both daytime and nighttime in Summer due to their similar wave characteristics and propagating directions.