1 Introduction

In the time before or after an earthquake, certain unusual natural phenomena can be observed, such as geothermic increases near the earthquake epicenter, abnormal increases of chemical gases (\(\hbox {O}_{3}\), \(\hbox {N}^{+}\), Rn, \(\hbox {CO}_{2}\) etc.), and evident ionospheric electron density disturbances (Cervone et al. 2006; Oyama et al. 2008; Pulinets and Boyarchuk 2004, Chapter 1.2.2–1.2.4). Among those anomalies, Co-Seismic Ionospheric Disturbances (CID) and Pre-Earthquake Ionosphere Abnormal (PEIA) have gradually attracted the attention of many scholars due to their relatively short-term manifestation immediately after the main shock or within 15 days before the occurrence of earthquake. The ionospheric anomalies associated with earthquakes first received attention in 1965 (Bolt 1964; Davies and Baker 1965; Leonard and Barnes 1965). Over the last few decades, many statistical analyses of seismic ionosphere anomalies were conducted to better understand this behavior (Liu et al. 2006; Kamogawa 2006; Saroso et al. 2008; Kon et al. 2011; Le et al. 2011; Yao et al. 2012; Perevalova et al. 2014). The electric field penetration of the lithosphere–atmosphere–ionosphere coupling (LAIC) is one of the most accepted interpretations for the Pre-Earthquake Ionosphere Abnormal (Pulinets and Davidenko 2014; Sorokin et al. 2007; Freund et al. 2009; Ryu et al. 2014a). A possible mechanism is noted by Sorokin et al. (2007). The release of soil gases accompanied with rock rupturing transfers positive and negatively charged aerosols into the atmosphere. Atmospheric convection of the charged aerosols forms external electric currents, which work as a source of perturbation in the atmosphere–ionosphere electric circuit. The abnormal electric field would cause the disturbance of electron density in the ionosphere due to \(E \times B\) force. This process is called LAIC electric field penetration. The recent simulation results from the SAMI3 model show that the applied zonal electric field leads to TEC variations near the source region, a westward electric field causes a downward \(E \times B\) motion of the ionosphere, and an eastward electric field causes an upward \(E \times B\) motion of the ionosphere (Kuo et al. 2014, 2015). Pulinets and Davidenko (2014) indicate that the acoustic gravity wave (AGW) mechanism has not shown its effectiveness for recent large earthquakes. But to date, it should be noted that the physical mechanisms of the PEIA remain unclear.

The AGW excited by earthquakes can propagate to ionospheric heights and generate traveling ionospheric disturbances (TIDs). The AGW is considered to be a possible physical mechanism that triggers CID (Liu et al. 2010; Heki 2011; Maruyama et al. 2011). In 1995, Calais and Minster (1995) noted that CID are caused by the propagation of AGW when they studied the Northbridge earthquake ionospheric disturbance. This is the first case in which a CID was investigated by monitoring the total electron content (TEC). AGW-ionosphere coupling was first proposed by Yuen et al. (1969) with ionosonde data. Afraimovich et al. (2010) also found that the wave propagation direction of ionospheric disturbances is consistent with the earthquake shattering direction by studying the ionospheric disturbance of the Wenchuan earthquake. The subject of earthquake ionosphere anomalies is a discipline filled with challenges and could support the classic techniques for earthquake and tsunami warning system (Jin et al. 2015).

At UTC 06:11:26.9 (LT11:56:26.9) on April 25, 2015 (Day of Year, DOY115), a \(\hbox {Mag}=7.8\) earthquake occurred in the Nepal region, with an epicenter located at 28.24\(^{\circ }\)N, 84.74\(^{\circ }\)E and a focal depth of 15 km. Reddy and Seemala (2015) reported the ionospheric response induced by the acoustic waves and Rayleigh waves of this earthquake and noted that the disturbances were distributed among certain areas (zone 1). Kamogawa et al. (2015) also found an ionospheric hole caused by this large inland earthquake and argued that an ionospheric hole could also be observed after a significantly large inland earthquake, not just following submarine earthquakes. Sun et al. (2016) found that the vertical scale of the earthquake ionosphere disturbance is 150 km, while the \(\hbox {h}_{\mathrm{m}}\hbox {F}_{2}\) is uplifted by more than 30 km within 1 min. Chen et al. (2016) reported the ionosphere disturbances of this earthquake as far as 1500–2500 km to the northeast of the epicenter. This paper mainly focuses on the Co-Seismic Ionospheric Disturbances of the Nepal earthquake around the main shock time and investigates the potential physical processes by data analysis and simulation results. In Sect. 2, we present the geomagnetic conditions (Dst index, Kp index) and solar radiation EUV (F10.7) for 30 days before and after the earthquake, as the fluctuations of these factors could cause major ionospheric disturbances. In particular, we analyze the geomagnetic conditions and solar activity around the main shock time. In Sect. 3, the TEC estimation method is introduced, and the estimation accuracy of differential code bias (DCB) product is given to demonstrate the reliability of the TEC sequence used in this paper. Then, in Sect. 4, by using several arcs of the Satellite-Station TEC (ssTEC) sequence, the occurrence time and spatial location of the CID are determined. Knowing the tridimensional reconstruction of the CID may improve the early warning system of earthquakes and tsunamis (Occhipinti et al. 2008), so in Sect. 5, the CID process of specific time periods and locations is tracked and analyzed by CIT images. In final discussion section, the potential physical mechanism is investigated by using the simulation experiment with seismic electric penetration, which confirms the CIT results.

2 Geomagnetic, solar, and meteorological conditions

The influence intensity of the geomagnetic conditions and solar radiation in the ionosphere at a certain time point can be reflected by three types of index variation sequences: the Dst index, the Kp index, and the F10.7 index. Figure 1 shows three index variation sequences for 30 days before and after the main shock time. First, on the earthquake day (DOY115), the Dst index variations are within − 10 nT, and the Kp indexes are below 1. This indicates that the geomagnetic conditions were extremely quiet on that day. According to the F10.7 index variations on DOY115, the solar EUV radiation is also in a degraded state, within 1.5% of its mean value. The three indexes taken together show there are no abnormal changes on DOY115.

In addition, it can be seen from DOY101 to DOY114 that the Dst index experiences two troughs (at DOY101 and DOY106), with orders of magnitude of approximately − 100 nT, and the Kp index shows four peaks, with all peak values beyond 4. During this period, the F10.7 index exhibits an increasing trend at the beginning and then remains at the maximum level with an order of magnitude above 150 until DOY111. After that, the F10.7 index begins to decrease.

In brief, the geomagnetic and solar conditions are relatively quiet and stable from DOY91 to DOY100. The Dst index remains within − 30 nT, and the Kp index remains less than 3, except for a maximum value on DOY92. According to the determination standard value for a geomagnetic storm, no geomagnetic storm occurred during this period. The F10.7 index also remains at a low magnitude. The Dst and Kp indexes on DOY114 and DOY116 are also at lower levels. Thus, the observations on these 12 days are used to derive the differential code bias (DCB) products and background values.

According to the meteorological data recorded by the Station Kathmandu Airport (27.70\(^{\circ }\)N, 85.37\(^{\circ }\)E), there was no rain or melted snow before the earthquake from DOY91 to DOY115, except on DOY94 (0.12 in), DOY95 (0.13 in), DOY102 (0.01 in), DOY103 (0.39 in), DOY108 (0.4 in), and DOY112 (0.13 in). In particular, on the two days before, there was no rain around this area. Based on that, the ionosphere disturbances could not be caused by the extreme weather conditions.

3 The dataset and estimation method of TEC and DCB

Due to the dispersion effect, the GNSS dual-frequency observation data are used to estimate the ionosphere electron density. Because TEC can be estimated with high accuracy and resolution, GNSS ionospheric inversion has become an important means for monitoring and studying the ionosphere. Formulas (1)–(4) provide the basic GNSS observation equations:

$$\begin{aligned} L_1= & {} \lambda _1 \Phi _1 \nonumber \\= & {} \rho -I+\lambda _1 N_1 +b_1^T +b_1^R +M(\Phi _1 )+\varepsilon (\Phi _1 ) \end{aligned}$$
(1)
$$\begin{aligned} L_2= & {} \lambda _2 \Phi _2 \nonumber \\= & {} \rho -\gamma I+\lambda _2 N_2 +b_2^T +b_2^R +M(\Phi _2 )+\varepsilon (\Phi _2 ) \end{aligned}$$
(2)
$$\begin{aligned} P_1= & {} \rho -I+B_1^T +B_1^R +M(P_1 )+\varepsilon (P_1 ) \end{aligned}$$
(3)
$$\begin{aligned} P_2= & {} \rho -\gamma I+B_2^T +B_2^R +M(P_2 )+\varepsilon (P_2 ) \end{aligned}$$
(4)

where \(\Phi _{1} \), \(\Phi _{2} \), \(P_{1} \), and \(P_{2} \) are the two carrier phase observations and two pseudo-range observations, respectively, and \(L_{1} \) and \(L_{2} \) are the carrier observations with units of meters. \(\lambda _{1} \) and \(\lambda _{2} \) are the wavelengths of \(L_{1} \) and \(L_{2} \), respectively, and \(\rho \) is the geometric distance between the receiver and satellite, including the errors independent on frequency. I is ionospheric delay, and \(\gamma \) is the ratio parameter, expressing the relationship between the ionospheric delays of the two carrier observations. \(\gamma ={f_1^2 }/{f_2^2 },(1-\gamma )=-0.6469\). \(N_1 \) and \(N_2 \) are ambiguity parameters. \(B_k^T \) and \(B_k^R \) are satellite and the receiver code biases, and \(b_k^T \) and \(b_k^R \) are satellite and the receiver carrier phase biases, respectively. T and R represent the satellite and receiver, respectively, and \(k=1,2\) represents two different frequencies. M and \(\varepsilon \) are multipath error and measurement noise, respectively.

The difference between the two frequency observation equations [i.e., (1) and (2), and (3) and (4), respectively] is as follows:

$$\begin{aligned} L_1 -L_2= & {} \left( {1-\gamma } \right) I_1 +(\lambda _1 N_1 -\lambda _2 N_2 )+\left( {b_1^T -b_2^T } \right) \nonumber \\&+\left( {b_1^R -b_2^R } \right) +\left( {M_{L1} -M_{L2} } \right) +\varepsilon _{L1,2} \end{aligned}$$
(5)
$$\begin{aligned} P_1 -P_2= & {} \left( {1-\gamma } \right) I_1 +\left( {B_1^T -B_2^T } \right) +\left( {B_1^R -B_2^R } \right) \nonumber \\&+\left( {M_{p1} -M_{p2} } \right) +\varepsilon _{p1,2} \end{aligned}$$
(6)

where \(\varepsilon _{1,2} \) is the difference of two frequency noise levels. \(B^{T}\) and \(B^{R}\) are differential code biases (DCBs), where \(B^{T}=B_1^T -B_2^T \), \(b^{T}=b_1^T -b_2^T \), \(B^{R}=B_1^R -B_2^R \), and \(b^{R}=b_1^R -b_2^R \).

To reduce the impact of multipath effects on TEC estimation, the appropriate cutoff angle (\(15^{\circ }\)) is selected. For dual-frequency GPS observations, the ionospheric delay of the L1 frequency can be expressed as:

$$\begin{aligned} I_1 =\frac{40.3}{f_1^2 }\hbox {TEC} \end{aligned}$$
(7)

Substituting (7) into (5) and (6), the TEC can be expressed as:

$$\begin{aligned} \hbox {TEC}_p= & {} \frac{f_1^2 f_2^2 }{40.3\left( {f_2^2 -f_1^2 } \right) }\left[ {\left( {P_1 -P_2 } \right) -B^{T}-B^{R}} \right] \end{aligned}$$
(8)
$$\begin{aligned} \hbox {TEC}_L= & {} \frac{f_1^2 f_2^2 }{40.3\left( {f_2^2 -f_1^2 } \right) }\nonumber \\&\left[ {\left( {L_1 -L_2 } \right) -(\lambda _1 N_1 -\lambda _2 N_2 )-b^{T}-b^{R}} \right] \end{aligned}$$
(9)

It can be observed from Formulas (8) and (9) that reliable DCB estimations are the prerequisite to guaranteeing high-accuracy TEC. To more accurately extract CID fluctuations on the earthquake day (DOY115), when computing the TEC sequence on DOY115, the fixed satellite and receiver DCBs are adopted in the estimation process. According to Fig. 1, the observation data from DOY91 to DOY100 are used for the DCB estimation in advance, as the geomagnetic conditions were quiet on those days. Considering that the DCB parameters are stable over a month, the average DCB values of 12 days are used in the TEC estimation on DOY115.

As shown in Fig. 2, GNSS data used in this paper are from 12 stations of the Crustal Movement Observation Network of China in Tibet and 10 stations of the Nepal GNSS network. The data sampling interval is 30 s.

Fig. 1
figure 1

Dst, Kp, F10.7 index and rains record variations from DOY91 to DOY121 of 2015 (the red line marks the earthquake day)

Full size image
Fig. 2
figure 2

The distribution of the stations used (red star marks the epicenter)

Full size image

It can be seen from the table that the standard deviations are mostly within 0.5 ns. The table gives the average values for each receiver, and the values will be used for the subsequent TEC estimation as fixed values (Table 1).

Table 1 Average values and standard deviation of the 22 receivers’ DCB parameters (unit: ns)
Full size table

4 The principle and reliability 1D CID detection method

The measurement accuracy of carrier observation is 100 times higher than the code pseudo-range, so the carrier observation is used here to analyze the TEC disturbances. Processing the ssTEC sequence with the second-order operator is an effective method for CID detection, as the ambiguity parameters are subtracted with the observation differences between epochs. The principle of this method is to predict the current epoch’s normal ionospheric changes using the prior two epochs’ TEC variation. After removing the normal ionospheric changes of the current epoch, the abnormal changes can be detected precisely. The cycle slip is detected with the TurboEdit method in advance.

Fig. 3
figure 3

The ssTEC sequences of Satellite 23 observed from different stations around the seismic time (scale in the figure is 0.3 TECU)

Full size image

As shown in Formulas (10) and (11), here, we have a set of ssTEC time sequences:

$$\begin{aligned} \hbox {TEC}_1 ,\ldots ,\hbox {TEC}_{i-1} ,\hbox {TEC}_i ,\hbox {TEC}_{i+1} ,\ldots ,\hbox {TEC}_n \end{aligned}$$
(10)

The ionospheric disturbance of the ith epoch could be determined by the TEST value given by Formula (11):

$$\begin{aligned}&\hbox {TEC\,TEST}=\hbox {TECT}=\hbox {TEC}_i -(\hbox {TEC}_{i-1} +\hbox {TECV}_{\mathrm{EST}} )\nonumber \\&\hbox {TECV}_{\mathrm{EST}} =\hbox {TEC}_{i-1} -\hbox {TEC}_{i-2}=\frac{f_1^2 f_2^2 }{40.3(f_2^2 -f_1^2 )} \nonumber \\&\quad [\Delta L_1 (i-1)-\Delta L_2 (i-1)] (\mathrm{see~formula}~9) \nonumber \\&\mathrm{Then} \nonumber \\&\hbox {TECT}=\frac{f_1^2 f_2^2 }{40.3(f_2^2 -f_1^2 )}[(\Delta L_1 (i)-\Delta L_2 (i))\nonumber \\&\qquad \qquad -(\Delta L_1 (i-1)-\Delta L_2 (i-1))] \nonumber \\&\Delta L_j (i)=L_j (i)-L_j (i-1),\quad j=1,2 \end{aligned}$$
(11)

where TECT (TEC TEST) is the test value for the disturbances detection. \(\hbox {TECV}_{EST} \) is the estimated TEC change velocity between the adjacent epochs. Combined with the ionosphere pierce point (IPP) coordinates of the TEC sequence, the 1D CID detection method can provide the CID time and location of each observation arc. The occurrence time and location are useful references for the next 2D or 3D detection. Here the ionosphere altitude for single-layer model is set at 450km.

Figure 3 shows the ssTEC sequences of Satellite 23 observed from different stations. The red line in the figure marks the main shock time. From the figure, a clear ionospheric disturbance can be observed among the seven ssTEC sequences, corresponding to BRN2, RMTE, TPLJ XZAR, XZRK, XZZB, and XZZF, approximately 10 min after the main shock. Among the stations, the reaction amplitudes of BRN2, RMTE, and TPLJ are much smaller than those of the other four stations. According to the IPP distribution given in Fig. 4, however, the IPP coordinates of BRN2, RMTE, and TPLJ are closer to the epicenter than those of the other four stations.

The ssTEC sequences of XZCY and XZCD experience the TEC fluctuations during the period of interest. Additionally, there are no obvious disturbances among other stations’ ssTEC sequences after the main shock. This phenomenon is consistent with the findings of Reddy and Seemala (2015).

Fig. 4
figure 4

IPP distribution of Satellite 23 for the 22 stations (the black \(\times \) marks the IPP locations at the time when earthquake occurred, and the blue \(\times \) marks the IPP locations where the TEC sequences show obvious ionospheric disturbances)

Full size image
Fig. 5
figure 5

The TEC variation of Satellite 3 of different stations around the seismic time (scale in the figure is 0.3 TECU)

Full size image

To further confirm the CID reaction area, in Fig. 5, the ssTEC sequences of Satellite 3 are given. Compared with Fig. 3, among the seven stations presenting the CID, stations RMTE and XZZB fail to show the CID reaction in the Satellite 3 ssTEC sequence. According to the IPP coordinates in Fig. 6, the IPPs shift to the west of the epicenter, especially to station XZZB; the IPPs of Satellite 3 are far away from the epicenter compared with those of Satellite 23. To summarize, according to the ssTEC sequences, the influences of CID are shown only in the stations in the special area (zone 1 in Reddy and Seemala 2015) near the epicenter, but the disturbances’ amplitude is not proportional to the distances between the stations and the epicenter. Also in Fig. 5, as Kamogawa et al. (2015) stated in their paper, the ionospheric hole could be clearly distinguished in sequences of BRN2 and TPLJ.

Fig. 6
figure 6

IPP distribution of Satellite 3 for 22 stations (the black \(\times \) marks the IPP locations at the time when the earthquake occurred, and the blue \(\times \) marks the IPP locations where TEC sequences show obvious ionospheric disturbances)

Full size image

As the measurement accuracy of carrier phase observations is at the millimeter level, the timing and location of the ionospheric disturbance given by the 1D CID detection is credible. As shown in Formulas (8) and (9), compared with the pseudo-range observations, the carrier phase has the problem of the uncertainty of the ambiguity. In Formula (11), Formula (9) considers the difference between epochs, and the ambiguity can be removed if no cycle slip occurs. Hence, using formula (11), the TEC variation between epochs can be estimated. As shown in Formula (12):

$$\begin{aligned} \hbox {TECV}_{\mathrm{EST}} =\frac{f_1^2 f_2^2 }{40.3\left( {f_2^2 -f_1^2 } \right) }\left[ {\Delta L_1 -\Delta L_2 } \right] \end{aligned}$$
(12)

Therefore, the estimation precision of \(\hbox {TECV}_{\mathrm{EST}} \) between epochs can be obtained from Formula (13):

$$\begin{aligned}&\sigma _{\mathrm{TECV}_{\mathrm{EST}} }^2 =\left[ \frac{f_1^2 f_2^2 }{40.3\left( {f_2^2 -f_1^2 } \right) }\right] ^{2}\left[ {\sigma _{\Delta L_1 }^2 +\sigma _{\Delta L_2 }^2 } \right] ,\nonumber \\&\frac{f_1^2 f_2^2 }{40.3\left( {f_2^2 -f_1^2 } \right) }=-{9}.{5196} \nonumber \\&\sigma _{L1} =\lambda _1 /100=0.0019m;\sigma _{L2} =\lambda _2 /100=0.0024m \nonumber \\&\sigma _{\mathrm{TECV}_{\mathrm{EST}} }^2 =(-{9}.{5196)}^{2}(2\times 0.0019^{2}+2\times 0.0024^{2}) \nonumber \\&\sigma _{\mathrm{TECV}_{\mathrm{EST}} } =0.04\,(\hbox {TECU}) \end{aligned}$$
(13)

where \(\sigma \) is the standard deviation. From Formula (13), it can be observed that the \(\hbox {TECV}_{\mathrm{EST}} \) accuracy is approximately 0.04TECU. Here we must illustrate that the CIT results in Fig. 7 are average images over 10 min, and the TEC disturbances in Figs. 3 and 5 show TEC sequences between two neighbor epochs. Additionally, the background values are different in the two methods: in Figs. 3 and 5, the background is the prediction value based on the two previous observations, while in Fig. 7, the background is the average TEC values on the same time span of DOY114 and 116. For these two reasons, the disturbance amplitudes of the figures are different, but the disturbance occurrence times of the two methods are consistent.

5 CIT images

To extract the fluctuation information, the background values should be selected first. According to Fig. 1, the geomagnetic conditions and solar radiation were quite calm on DOY114 and DOY116. Here, the cutoff angle is set at 15\(^{\mathrm{o}}\), and the average values of TEC on DOY114 and DOY116 are taken as the background values. The percentage increase can be calculated with Formula (14),

$$\begin{aligned}&\hbox {TEC}_{114} ,\hbox {TEC}_{115} ,\hbox {TEC}_{116} \nonumber \\&\hbox {TEC}_{\mathrm{aver}} =\frac{(\hbox {TEC}_{114} +\hbox {TEC}_{116} )}{2} \nonumber \\&\% \hbox {dif}=\frac{(\mathrm{TEC}_{115} -\mathrm{TEC}_\mathrm{aver} )}{\mathrm{TEC}_\mathrm{aver} } \end{aligned}$$
(14)
Fig. 7
figure 7

CIT images around the seismic area and the main shock time

Full size image

GNSS-based computerized ionospheric tomography (CIT) is an effective method to study the 3D structure of the seismic ionospheric disturbances. As CIT is a typical ill-posed problem, it is hard to obtain reliable 3D ionosphere variation. There are two main methodological categories for solving this problem: the pixel-based method (PBM) (Avinash and Malcolm 1988; Pryse and Kersley 1992) and the function-based method (FBM) (Hansen 1998; Schmidt et al. 2008). Compared with PBM, FBM can acquire the continuously distributed electron density, while the solution is more unstable. Another problem is that the resolution time of CIT algorithms is too large to capture the detailed CID variation. Kong et al. (2016) introduced a so-called MFCIT (CIT with mapping function) algorithm. MFCIT splits the integration process into 4 layers (150, 450, 750, and 1050 km) along the observation ray, and then, the single-layer model (SLM) is applied to each integration part using a mapping function. The model parameters are estimated layer by layer with the Kalman filtering method. In this way, the resolution time can be decreased to 10 min and the CIT accuracy could be ensured but with the space resolution decreased. So the new algorithm allows us to capture the large-scale vertical change of ionosphere during the disturbances. From Figs. 3 and 5, it is seen that the disturbances measured by some stations lasted for over 10 min, so in this paper, the MFCIT is used to invert the 3D electron density distribution. After the TEC distributions are estimated using MFCIT, the TEC increase percentages at different layers are computed using Formula (14). Figure 7 shows an electron density increase percentage on the 150- and 450-km layers from UT06:00 to UT06:30 with 10-min resolution, at a local time of LT11:45 to LT12:15.

In Fig. 7, from top to bottom, the CIT time periods are UT06:00-UT06:10, UT06:10-UT06:20, and UT06:20-06:30, respectively. Based on the right three images, the TEC variation on the 450-km layer shows no anomaly before or after the main shock, while the three images under different time frames present a similar structure. On the northeast of the images, there is a decrease hole, which may be related to the topographic factors. As the edge of Qinghai–Tibet Plateau ends here, the great altitude difference may trigger tangential wind and then cause ionospheric disturbances. It could be seen from the TEC distribution on the 150-km layer under the same period the similar decreasing anomaly can be found. To the southwest, there is an increasing anomaly, which may be related to solar activity; according to the local time period (LT11:45 to LT12:15), direct sunlight happened to pass this area. Based on the left-hand images of the 150-km layer, a similar increasing anomaly can be found.

During the third time period (UT06:20-06:30), however, a difference is apparent between the TEC distribution at 150 km and 450 km. The abnormal TEC disturbances are shown in the third image in the 150-km layer (UT06:20-06:30). It can be seen that on the left of the epicenter, the north area shows increasing disturbances, while the south shows decreasing disturbances. On the right of the epicenter, the opposite phenomena can be found. According to Figs. 3 and 5, the 2D TEC disturbances are first observed at \(\sim \)UT0620 and last for \(\sim \)10 min. The 3D images also show the disturbances during this period, but the anomalies are limited to 150 km. In fact, the penetration altitude of the anomalous electric field is also related to the geomagnetic location, and the magnetic field intensity will influence the penetration altitude. The physical mechanism will be discussed in next section.

6 Discussion

Astafyeva et al. (2013) demonstrated that shallow earthquakes with magnitudes of 7.2 M–7.8 M could cause co-seismic perturbations with near-field amplitudes of 0.2–0.4 TECU (lasting 4–8 min), whereas mega earthquakes of \(\sim \)Mw 9.0 and above produce extremely large perturbations of \(\sim \)1–3 TECU (lasting 30–40 min). As no evident wave-like disturbances were presented by the GPS-TEC observation, we consequently assume that the seismic electric penetration could be a possible interpretation.

Fig. 8
figure 8

The distribution of the abnormal electric field (mV/m) around the epicenter in the ionospheric region (the red triangle represents the epicenter)

Full size image
Fig. 9
figure 9

The ionospheric electron density disturbances in the geomagnetic plane west of the epicenter and east of the epicenter, which is calculated by the SAMI2 (the red point represents the epicenter, and the latitude of the west plane and east plane is 81\(^{\mathrm{o}}\)E, 91\(^{\mathrm{o}}\)E, respectively)

Full size image

PEIA is generally interpreted as the lithosphere–atmosphere–ionosphere coupling (LAIC) mechanism (Pulinets et al. 2000; Sorokin et al. 2001). The most important LAIC process is the penetration of an abnormal electric field into the ionosphere, which consequently induces \(E \times B\) disturbances. Here we utilize a three-dimensional electric field penetration model for seismo-ionospheric disturbance (Zhou et al. 2017). This model is similar to the thunderstorm electric field penetration model, in which the electric potential changes in response to the abnormal electric field penetration are solved in the dipole coordinates. By considering the magnetic conjugation effects, the electric potential is integrated along the magnetic field lines from the lower boundary of the southern hemisphere to that of the northern hemisphere. The model is capable of producing an abnormal electric field in the ionosphere due to the LAI electric field penetration. Figure 8a presents the simulation results of the abnormal ionospheric electric field distribution in the XZ plane, where X is the direction of the magnetic field west to east and Z is the vertical direction to the ground. The parameters utilized in the simulation are adopted from the international reference ionosphere model. The white arrows show the intensity and direction of the penetrated electric field. Figure 8b presents the simulation results of the abnormal ionospheric electric field distribution in the XY plane, where Y is the direction of the magnetic field south to north. This simulation result indicates that an east to west abnormal electric field exists in the ionospheric region above the epicenter. The maximum electric field is approximately 4 mV/m and exists between 150 and 200 km, which can consequently cause notable redistribution of the ionospheric density due to the \(E\times B\) effects. Due to the magnetic field line effects, the penetration height of the abnormal electric field generated in the ionosphere decreases with the geomagnetic latitude. In this study, the geomagnetic latitude of the epicenter is 19.21\(^{\circ }\)N and the effective penetration height of the abnormal electric field is less than 400 km in our simulations, which is generally consistent with the CIT results reported in previous section.

We also reproduce the ionospheric electron density disturbances in the vicinity of the earthquake by using the SAMI2 model. SAMI2 calculates dynamic plasma and chemical evolution of ions and electrons in the altitude range from 100 km to several thousand kilometers in the middle and low latitudes (Huba et al. 2000). By integrating the abnormal electric field into the SAMI2 model, the re-distribution of the ionosphere electron density is presented in Fig. 9. From Fig. 9, it can be clearly observed that electron density increases in the north and decreases in the south in the west magnetic meridian plane of the epicenter, while the electron density decreases in the north and increases in the south in the east magnetic meridian plane of the epicenter. According to Fig. 8b, the direction of the abnormal electric field is the opposite, which generates different directions of electron movement.

7 Summary

In this paper, we use observational data from the China Mainland Construction Environmental Monitoring Network and Nepal Network to determine the timing and location of the CID. The methodology of TEC variation estimation between epochs obtained by the carrier phase for TEC anomaly detection is also introduced. According to the ssTEC sequences of Satellites 3 and 23 for 22 stations, the CID occurred approximately 10 min after the earthquake. The geomagnetic indexes (Dst, Kp index), solar activity parameters (F10.7 index), and the meteorological condition demonstrate the geomagnetic conditions, solar radiation, and other factors that can cause ionospheric disturbances were calm on that day. Meanwhile, only seven of 22 stations show obvious ionospheric disturbances, which means the ionospheric disturbance was a regional reaction. The ionospheric hole is also distinguished from the Satellite 3 ssTEC sequence of the two Nepal network stations.

Using data from 22 stations near the epicenter, we reconstruct three-dimensional ionosphere images with resolution of 10 min. Using the electron density on DOY114 and DOY116 as the background values, it is found that the seismic ionospheric disturbance of this Nepal earthquake was limited to the 150-km layer, while the 450-km layer failed to present the anomaly. The penetration height of LAI electric field in the ionosphere is lower at low latitudes than at high latitudes, with the simulation experiment of electric field penetration in area of interest. It is found that the effective penetration height of the abnormal electric field is less than 400 km, which confirms the CIT results.

Additionally, we found that on the western side of the epicenter, the electron density was enhanced in the north but weakened in the south. On the eastern side of the epicenter, the opposite phenomenon could be observed, i.e., electron density was enhanced in the south but weakened in the north. Our main findings are summarized as follows.

  1. (1)

    The ssTEC results indicate that an evident depletion of ionospheric electron density is observed at the epicenter \(\sim \) 10 min after the earthquake.

  2. (2)

    Three-dimensional CIT results reveal that a) the ionospheric electron density disturbance primarily occurred below the 400 km height, and b) the west side and east side of the epicenter presented opposite distributions of the ionospheric electron density.

  3. (3)

    We introduce an electric field penetration model for seismo-ionospheric effects whose the abnormal electric field in the ionosphere is integrated to the SAMI2 model. The simulation results are consistent with the observations and generally interpret the potential physical process.