The very low-frequency transmitter radio wave anomalies related to the 2010 Ms 7.1 Yushu earthquake observed by the DEMETER satellite and the possible mechanism

Earthquakes may disturb the lower ionosphere through various coupling mechanisms during the seismogenic and coseismic periods. The VLF (very low-frequency) signal radiated from ground-based transmitters will be affected when it penetrates the disturbed ionosphere above the epicenter area, and this anomaly can be recorded by low-Earth orbit satellites under certain conditions. In this paper, the temporal and spatial variation of the signal-to-noise ratio (SNR) of the VLF transmitter signal in the ionosphere over the epicenter of 2010 Yushu Ms 7.1 earthquake in China is analyzed using DEMETER (Detection of Electro-Magnetic Emission Transmitted from Earthquake Regions) satellite observation. The results show that SNR over the epicenter of the Yushu earthquake especially in the southwestern region decreased (or dropped) before the main shock, and a GPS–TEC (Global Positioning System; total electron content) anomaly accompanied, which implies that the decrease in SNR might be caused by the enhancement of TEC. A full-wave method is used to study the mechanism of the change in SNR before the earthquake. The simulated results show SNR does not always decrease before an earthquake. When the electron density in the lower ionosphere increases by 3 times, the electric field will decrease about 2 dB, indicating that the disturbed-electric-field decrease of 20 % compared with the original electric field and vice versa. It can be concluded that the variation of electron density before earthquakes may be one of the important factors influencing the variation of SNR.

1 Introduction

The VLF (very low-frequency) radio waves radiated by the powerful ground-based VLF transmitters have been used for long-distance communication and submarine navigation because of the efficient reflection within the Earth–ionosphere waveguide. However, there is still a small fraction of the wave energy that can leak into the higher ionosphere and magnetosphere after being absorbed intensively by the lower ionosphere. The signals from transmitters observed by the LEO (low-Earth orbit) satellites can be used to study the propagation of VLF waves in the Earth–ionosphere waveguide and ionosphere, as well as wave–particle interaction in the radiation belt (Inan et al., 2007; Inan and Helliwell, 1982; Lehtinen and Inan, 2009; Parrot et al., 2007).

It is gradually confirmed that earthquake precursors not only appear near the ground but also may couple with the atmosphere and ionosphere through some mechanisms, resulting in plasma disturbances in the ionosphere and recorded by various instruments like ionosonde or GPS (Global Positioning System) receivers measuring TEC (total electron content) (Liu et al., 2009, 2001, 2006; Pulinets et al., 2000; Stangl et al., 2011; Zhao et al., 2008). Therefore, the amplitude of the VLF signals from the ground-based VLF transmitter observed on the ground and from satellites will change when encounter the disturbed area in the ionosphere (Hayakawa, 2007; Maurya et al., 2016; Molchanov et al., 2006; Píša et al., 2013). Molchanov et al. (2006) have found the signal-to-noise ratio (SNR) of the electric field from VLF transmitters recorded by the DEMETER (Detection of Electro-Magnetic Emission Transmitted from Earthquake Regions) satellite decreased near the epicenters during a series of earthquakes. The spatial size of an SNR reduction zone increases with the magnitude of the earthquake. However, it is hard to distinguish the coseismic anomaly and precursor from their results.

Two devastating earthquakes, the 2008 Ms 8.0 Wenchuan earthquake and the 2010 Ms 7.1 Yushu earthquake, occurred successively in southwestern China during the operation period (2004–2010) of the DEMETER satellite. Some research have also focused on the SNR variation of VLF transmitters using DEMETER satellite observation to extract the earthquake-related anomalies before the two strong earthquakes (He et al., 2009; Shen et al., 2017; Yao et al., 2013). The results all illustrated the decrease in SNR before the earthquakes. Since the earthquake-related ionospheric disturbance zone is not right over the epicenter, the relative position of the SNR anomaly and the epicenter should be further studied. The factors which influence SNR and the possible mechanism also needs to be comprehensively illustrated.

The Alpha VLF transmitters in Russia transmit three frequencies in each station which provide us opportunities to study the influence of the ionosphere on different wave frequencies. The devastating earthquake nearest the transmitters in China is the 2010 Ms 7.1 Yushu earthquake. In this paper we investigate the temporal and spatial SNR variation of the VLF transmitter signal in the ionosphere near the epicenter of the Yushu earthquake using DEMETER observation. The background variations of SNR in the same period of 2007–2010 have also been studied to distinguish whether the SNR reduction is caused by an earthquake or just ionospheric background changes. The mechanism of how the seismo-ionospheric disturbance affects the variation of SNR is discussed in this paper.

Regarding the mechanism of the VLF radio wave variations in the altitude of a LEO satellite (presented as SNR variation) before the earthquakes, Hayakawa (2007) and Píša et al. (2013) suggest the VLF anomalies exist because the lower ionosphere is lowered before earthquake. Molchanov et al. (2006) declared that the variation of SNR of satellite data is attributed to the ionospheric disturbance, especially the lower-ionospheric disturbance. Furthermore, it has been found that the electron density variation could exist in the lower ionosphere according to the computer ionosphere tomography (CIT) results based on GPS–TEC data before the Nepal Ms 8.1 earthquake in 2015 (Kong et al., 2018). The electric-field-penetrating model has shown that the electron density and height of the lower ionosphere can be changed by the additional current in the global electric circuit before the earthquake. On the other hand, Marshall et al. (2010) construct a 3D finite-difference time domain model to simulate that lightning could also cause the disturbance of the electron density in the lower ionosphere, which has a similar mechanism to the earthquake. Many studies also have found that the main loss of VLF wave power mainly occurs in the D–E region of the ionosphere when the wave penetrates into the ionosphere (Cohen and Inan, 2012; Liao et al., 2017; Starks et al., 2008; Tao et al., 2010; Zhao et al., 2017, 2015). In sum, the electron density variation in the lower ionosphere might be one main factor causing the SNR anomaly of VLF transmitter signal in the ionosphere. Based on these results, the full-wave calculation model was utilized to study the influence of the electron density disturbance of the lower ionosphere on the variation of VLF radio signals.

In this paper, a brief description of the DEMETER data and full-wave method used in this study are presented in Sect. 2. The temporal and spatial variations of SNR over the epicenter have been investigated before the Yushu earthquake with 4 years (2007–2010) of data; the full-wave model is used to simulate how the variation of electron density in the lower ionosphere affects the SNR of the electric field from a VLF transmitter at the altitude of a satellite in Sect. 3. The discussion and conclusions of this study are presented in Sects. 4 and 5 separately.

2 Materials and methods

2.1 Earthquake, VLF transmitters, and DEMETER data

At the local time of 07:49:37.9 on 14 April 2010, a Ms 7.1 earthquake hit the city of Yushu, Qinghai province, with an epicenter at 33.2^∘ N, 96.6^∘ E and a 14 km depth in the northeastern Tibetan Plateau. The nearest VLF transmitter around the epicenter is in the proximity of Novosibirsk (NOV), which belongs to the Russian Alpha navigation system consisting of three transmitters. The other two transmitters named Krasnodar (KRA) and Khabarovsk (KHA) are far away from the Yushu earthquake, so only the satellite data radiated from NOV have been used for analysis in this paper. The location of the transmitters and the epicenter of the Yushu earthquake are denoted by blue squares and black stars, respectively, in Fig. 1. Each transmitter radiates three different frequency VLF radio signals (11.9, 12.6, and 14.9 kHz), with a 0.4 s duration and a 3.6 s cycle.

Figure 1 The locations of transmitters and the Yushu earthquake. The blue squares represent the locations of the three transmitters (KRA, NOV, and KHA) in Russia. The epicenter of Yushu earthquake is denoted by the black star. The black square covers the region of the epicenter $\pm \mathrm{10}{}^{\circ }$, in which the data have been studied.

The DEMETER satellite was launched on 29 June 2004 as a sun-synchronous orbit at an altitude of 710 km, which then was changed to 660 km in December 2005 (Parrot et al., 2006), and the operation was ended in December 2010. The scientific objective of DEMETER is to detect and characterize the electromagnetic signals associated with natural phenomena (such as earthquakes, volcanic eruptions, and tsunamis) or anthropogenic activities. It operated in the region from invariant latitude −65 to 65^∘, with descending and ascending orbits crossing the Equator at the local time of ∼10:00 and ∼22:00, respectively. DEMETER has a revisit orbit period of about 14 d, which means the satellite returns over the same orbit trajectory after 13 d. The payloads include several electromagnetic sensors with two working mode: burst and survey. At the ELF–VLF (extremely low-frequency) band, the intensive-electromagnetic-wave data over locations of particular interest were provided in the burst mode, and in the survey mode, electric and magnetic power spectral density (PSD) data every 2 s were provided with a sampling frequency of 40 kHz and spectral resolution of 19.53 Hz.

According to the formula of Dobrovolsky et al. (1979) the preparation zone of the earthquake can reach ρ=10^0.43M, where M is the magnitude of the earthquake and ρ is measured in kilometers. Considering the limited extension of the Ms 7.1 Yushu earthquake, the preparation zone ρ can reach to 1130 km; we mainly focused on the region within the region of the epicenter $\pm \mathrm{10}{}^{\circ }$ (black square in Fig. 1). In this study, the nighttime PSD data of the electric field from the DEMETER's survey mode observations were extracted study the perturbations of the VLF signal before and after the Yushu earthquake. As the VLF radio signals at daytime are too small to cause obvious SNR variation compared with those at nighttime, we did not use the daytime data in this study.

2.2 The method to calculate SNR

According to the method of Molchanov et al. (2006), the SNR of the electric field was calculated as follows:

$$\begin{array}{}\text{(1)}& \text{SNR}={\displaystyle \frac{\mathrm{2}A\left(f_{\mathrm{0}}\right)}{A\left(f_{+}\right)+A\left(f_{-}\right)}},\end{array}$$

where A(f₀) is the amplitude of the electric-field spectrum at the central frequency and A(f_±) values are the spectrums at $f_{\pm }=f_{\mathrm{0}}\pm \mathrm{\Delta }f$, where Δf is the chosen frequency band. For the three Russian VLF transmitters, the f₀ is set as three VLF radio waves frequency radiated from NOV transmitters at 11.9, 12.6, and 14.9 kHz and the Δf=300 Hz.

2.3 Full-wave method

A full-wave method has been used to seek a solution of Maxwell equations for waves varying as e^jωt in a horizontally stratified medium with fixed dielectric permittivity tensors $\widehat{\mathit{\epsilon }}$ and permeability μ in each layer. Considering the region of our interest is much smaller than the radius of the earth, the earth's curvature is neglected in this study. A Cartesian coordinate system is established with x, y in the horizontal plane and z as vertical upward. We seek a solution of the Maxwell equations in the form of a linear combination of plane waves $\sim e^{j(k_{\mathrm{bot}}\cdot r_{\mathrm{bot}})}$, where k_bot is the horizontal component of the wave vector k, which is conserved by Snell's law inside each layer; we have

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}\mathit{k}\times E=\mathit{\omega }{\mathit{\mu }}_{\mathrm{0}}H\\ \mathit{k}\times H=-\mathit{\omega }\widehat{\mathit{\epsilon }}E\end{array}\right.,\end{array}$$

where ω is the angular frequency, μ is the permeability of the medium (µ≡1 for nonmagnetic medium), $\widehat{\mathit{\epsilon }}={\mathit{\epsilon }}_{\mathrm{0}}(I+\widehat{\mathit{\chi }})$ is dielectric tensor, and $\widehat{\mathit{\chi }}$ is electric-susceptibility tensor (Yeh and Liu, 1972). $\widehat{\mathit{\chi }}$ is determined by the electron density and collision frequency in the ionosphere, as well as the geomagnetic field. In our simulation, the electron density is calculated by the International Reference Ionosphere (IRI) model (Bilitza et al., 2017), and the electron collision frequency (denoted by v) is modeled by the exponential-decay law with the height (denoted by h) increasing at $v=\mathrm{1.8}\times {\mathrm{10}}^{\mathrm{11}}{\text{e}}^{-\mathrm{0.15}h}$. The parameters of the geomagnetic field at the location of the VLF transmitter is calculated by the International Geomagnetic Reference Field (IGRF) model (Finlay et al., 2010).

Eliminating the z components from Eq. (1), we can obtain the following elegant form of Maxwell equations.

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}z}}=j{\mathit{k}}_{\mathrm{0}}\widehat{T}\cdot V,\end{array}$$

where V=(E_botZ₀H_bot), Z₀ is wave impedance and $\widehat{T}$ is a 4×4 matrix:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& \widehat{T}=\\ & \left(\begin{array}{cc}\begin{array}{cc}-\frac{{\mathit{k}}_x{\mathit{\epsilon }}_{\mathrm{31}}}{{\mathit{k}}_{\mathrm{0}}{\mathit{\epsilon }}_{\mathrm{33}}}& -\frac{{\mathit{k}}_x{\mathit{\epsilon }}_{\mathrm{32}}}{{\mathit{k}}_{\mathrm{0}}{\mathit{\epsilon }}_{\mathrm{33}}}\\ -\frac{{\mathit{k}}_y{\mathit{\epsilon }}_{\mathrm{31}}}{{\mathit{k}}_{\mathrm{0}}{\mathit{\epsilon }}_{\mathrm{33}}}& -\frac{{\mathit{k}}_y{\mathit{\epsilon }}_{\mathrm{32}}}{{\mathit{k}}_{\mathrm{0}}{\mathit{\epsilon }}_{\mathrm{33}}}\end{array}& \begin{array}{cc}\frac{{\mathit{k}}_x{\mathit{k}}_y}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\epsilon }}_{\mathrm{33}}}& \mathrm{1}-\frac{{\mathit{k}}_x^{\mathrm{2}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\epsilon }}_{\mathrm{33}}}\\ -\mathrm{1}+\frac{{\mathit{k}}_y^{\mathrm{2}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\epsilon }}_{\mathrm{33}}}& -\frac{{\mathit{k}}_x{\mathit{k}}_y}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\epsilon }}_{\mathrm{33}}}\end{array}\\ \begin{array}{cc}-{\mathit{\epsilon }}_{\mathrm{21}}+\frac{{\mathit{\epsilon }}_{\mathrm{23}}{\mathit{\epsilon }}_{\mathrm{31}}}{{\mathit{\epsilon }}_{\mathrm{33}}}-\frac{{\mathit{k}}_x{\mathit{k}}_y}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{2}}}& -{\mathit{\epsilon }}_{\mathrm{22}}+\frac{{\mathit{\epsilon }}_{\mathrm{23}}{\mathit{\epsilon }}_{\mathrm{32}}}{{\mathit{\epsilon }}_{\mathrm{33}}}+\frac{{\mathit{k}}_x^{\mathrm{2}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{2}}}\\ {\mathit{\epsilon }}_{\mathrm{11}}-\frac{{\mathit{\epsilon }}_{\mathrm{13}}{\mathit{\epsilon }}_{\mathrm{31}}}{{\mathit{\epsilon }}_{\mathrm{33}}}-\frac{{\mathit{k}}_y^{\mathrm{2}}}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{2}}}& {\mathit{\epsilon }}_{\mathrm{12}}-\frac{{\mathit{\epsilon }}_{\mathrm{13}}{\mathit{\epsilon }}_{\mathrm{32}}}{{\mathit{\epsilon }}_{\mathrm{33}}}+\frac{{\mathit{k}}_x{\mathit{k}}_y}{{\mathit{k}}_{\mathrm{0}}^{\mathrm{2}}}\end{array}& \begin{array}{cc}-\frac{{\mathit{k}}_y{\mathit{\epsilon }}_{\mathrm{23}}}{{\mathit{k}}_{\mathrm{0}}{\mathit{\epsilon }}_{\mathrm{33}}}& \frac{{\mathit{k}}_x{\mathit{\epsilon }}_{\mathrm{23}}}{{\mathit{k}}_{\mathrm{0}}{\mathit{\epsilon }}_{\mathrm{33}}}\\ \frac{{\mathit{k}}_y{\mathit{\epsilon }}_{\mathrm{13}}}{{\mathit{k}}_{\mathrm{0}}{\mathit{\epsilon }}_{\mathrm{33}}}& -\frac{{\mathit{k}}_x{\mathit{\epsilon }}_{\mathrm{13}}}{{\mathit{k}}_{\mathrm{0}}{\mathit{\epsilon }}_{\mathrm{33}}}\end{array}\end{array}\right).\end{array}\end{array}$$

The electromagnetic field in each layer can be obtained in the k (wave vector) domain by solving Eq. (2) recursively in a direction which provides stability against numerical “swamping” (Budden, 1985; Lehtinen and Inan, 2008). The difficulty is how to deal with numerical stability when the solution of evanescent waves “swamp” the waves of interest because of the large imaginary vertical-wave number. More details of the full-wave method are described in Lehtinen and Inan (2008).

Figure 2 The evolution of SNR evolution VLF radio waves frequencies 11.9 kHz (top panel), 12.6 kHz (middle panel), 14.9 kHz (bottom panel) with Δf=300 Hz at nighttime. The black star stands for the epicenter of the Yushu earthquake; the gray line represents the transmitter that turns off on that day; the days with high geomagnetic activity are marked by blue and hollow dots. Please note that the date format of the black text in this figure is day-month-year, but the date format of the red and blue text is month day.

3 Results

3.1 VLF signal analysis from the DEMETER satellite

The SNR five revisit periods before and one revisit period after the earthquake in 2010 were calculated to study the evolution of SNR above the epicenter. The SNR distributions of three frequencies (11.9, 12.6, and 14.9 kHz) within the region of the epicenter $\pm \mathrm{10}{}^{\circ }$ are shown in Fig. 2, where the value of SNR is denoted with colored dots with different sizes and the black star represents the epicenter of the Yushu earthquake. The data during geomagnetic storms (here we defined Kp > 3 and Dst < −30 nT) were plotted with hollow dots, and very small gray dots mean that the transmitter is turned off on these days. It can be found that in the first revisit period (2–14 April) before the earthquake, the SNRs of the three frequencies all decrease dramatically compared with other periods no matter whether it is before or after the earthquake. In the first revisit period from 2 to 14 April in 2010, two magnetic storms occurred on 4–7 and 11–12 April, respectively.

Figure 3 The average SNR variation with revisit periods inside the square region with the center of the epicenter. Panels (a–c) show SNR at 11.9, 12.6, and 14.9 kHz and the numbers of the averaged data points. The green and red lines represent the SNR variations in 2010 and background time separately. The black dashed line represents the period with the end date of the main-shock date.

To minimize the impact of other factors and confirm whether the SNR anomaly is caused by the earthquake and not the variation of the ionospheric background, we focus on SNR in the black square (shown in Fig. 1) of the same period in 2007–2009 as the background, when there are no large earthquakes and the data when the transmitter was turned off or affected by geomagnetic storms are eliminated. The mean value of all the data in each period has been obtained to get the time sequence shown in Fig. 3. In Fig. 3, the black dashed line represents the occurrence date of the earthquake. The black and red solid lines represent the average values in the five periods before the earthquake and one period after the earthquake within the region of the epicenter $\pm \mathrm{10}{}^{\circ }$ in 2010 and background time, respectively. The change trends of SNR in the background time and 2010 are the same except in the first period before the earthquake. In the first period before the earthquake SNR decreased significantly in 2010, while it increased in background time at all transmitting frequencies. It means that the decrease in SNR in the first period in 2010 might be caused by Yushu earthquake.

The above results use the average value within the region of the epicenter $\pm \mathrm{10}{}^{\circ }$ in one revisit period of DEMETER to analyze the anomalies which ignore the day-to-day variability of the ionosphere. Furthermore, the daily variation of SNR in the first period before the earthquake is studied using a quartile-based process (Liu et al., 2009) to detect the anomaly of SNR. The median (M), the lower (first) quartile (denoted as LQ), and the upper (third) quartile (UQ) of every successive 11 d of SNR of the orbit data within the region of the epicenter $\pm \mathrm{10}{}^{\circ }$ have been calculated to find the deviation between the observed SNR of the 12th day and the computed median (M). Based on the assumption of the normal distribution of SNR with the mean (m) and standard deviation (σ), the expected value of M and LQ or UQ are equal to m and 1.34σ (Liu et al., 2009, and references therein). We set the lower boundary (LB in short) to LB = M + 2(M-LQ) and the upper boundary (UB) to UB = M + 2(UQ-M) to find the SNR anomalies with a stricter criterion. Thus, if an observed SNR on the 12th day is greater or smaller than its previous 11 d based on UB or LB, a positive or negative anomaly of SNR will be identified. Figure 4 shows the time series of SNR at 11.9, 12.6, and 14.9 kHz, and the red, gray, black curves denote the current SNR, associated median, and upper and lower boundary (UB and LB), respectively. Blue and green markers represent the positive and negative anomaly. As shown in Fig. 4, besides the negative anomalies which appeared on 13 April (1 d before the Yushu earthquake; the occurrence time of the Yushu earthquake is denoted by a vertical dashed line in Fig. 4) at all transmitting frequencies, another three anomalies occurred on 29 March, 8 April, and 10 April, respectively. Previous research indicates that earthquake anomalies usually occurred within 1 week before an earthquake, so the negative anomaly which occurred on 29 March at 12.6 and 14.9 kHz may be not related to the Yushu earthquake. The anomalies on 8 and 10 April only occurred on one single transmitting frequency, which may not be significant and is needed to be further researched.

Figure 4 A time series of SNR right above the Yushu epicenter. The Ms 7.1 Yushu earthquake occurred at the local time of 07:49:37.9 on 14 April 2010. The red, gray, and two black curves denote the currently observed SNR, associated median, and upper and lower bound (UB and LB), respectively. Blue and green sign represent the upper and lower anomalous days identified by the computer routine, respectively. The LB and UB are constructed by the previous 1–11 d moving median (M), lower quartile (LQ), and upper quartile (UQ), and the LB and UB are calculated by LB = M + 2(M-LQ) and UB = M + 2(UQ-M).

The result in Fig. 4 shows the anomalies of SNR during the successive 20 d before the Yushu earthquake. However, the 20 d orbital data may be carried into the ionospheric background noise of different space. To avoid this kind of ionospheric background noise, we select the three revisit orbits to analyze the anomalies of SNR before the Yushu earthquake further (the revisit orbit on 9 April overhead the epicenter, the revisit orbit on 13 April which is 550 km away from epicenter, and the revisit orbit on 10 April which is 750 km away from the epicenter are selected). The quartile-based process is also performed for every revisit's orbital data, but the 6 d sliding-mean value (including 3 d before the current day and 2 d after the current day) have been analyzed. The green and blue bar represent negative and positive anomalies in one orbit, respectively, in Fig. 5. As we can see in the top and middle panel, on 9 and 10 April, the negative and positive anomalies both occurred like other days in the same two revisit orbits. These anomalies could be induced by the daily variation. In the bottom panel, there are no obvious anomalies in other days with the same revisit orbit of 13 April, but the SNR values have obvious negative anomalies for all orbits on 13 April. These results further confirm that the anomalies of SNR occurred on 13 April.

Figure 5 A revisit orbital SNR of 9, 10, and 13 April 2010. The red, gray, and two black curves denote the currently observed SNR, associated median, and upper and lower bound (UB and LB), respectively. Blue and green bars represent the positive and negative anomalies in one orbit, respectively. The LB and UB are constructed by the 6 d moving median (M, including 3 d before the current day and 2 d after the current day), lower quartile (LQ), and upper quartile (UQ), and the LB and UB are calculated by LB = M + 2(M-LQ) and UB = M + 2(UQ-M).

We speculate that the anomalies of SNR may be related to the anomalies of electron density. To confirm our conjecture, we used GPS–TEC map data distributed by CODE (Center for Orbit Determination in Europe) to check out whether the total electron content (TEC) showed similar anomalies. The resolution of TEC data from CODE is $\mathrm{5}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$. We use 11 d sliding-mean values of every grid as a background, and then we can get a spatial distribution of the background. The background $\pm \mathrm{2}\times$ stand deviation is set as the threshold (upper bound and lower bound) to determine whether there are anomalies; if the intraday value exceeds the threshold, there are anomalies. We have reviewed the TEC anomalies of every day from 2 to 14 April (which means the duration of the sliding background is from 22 March to 13 April). The TEC anomalies only occurred on 13 April. The anomalies that are the most intensive were at 06:00 UT, which means only the SNR anomaly on 13 April is a possible earthquake precursor. The other two anomalies on 8 and 10 April in Fig. 4 may be caused by other factors. Figure 6a–b shows the TEC at 06:00 UT on 13 April and the sliding mean of the background (2–12 April); Fig. 6c–d shows the abnormal region where the TEC value exceeded the threshold (background $\pm \mathrm{2}\times$ standard deviation). As we can see, TEC had an abnormal enhancement on 13 April at the southwestern region of the epicenter. In addition, we collect the Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC) data in the abnormal region of TEC (southwestern region of the Yushu epicenter) to check whether there is abnormal variation in the D–E region electron density. As shown in Fig. 7, the result shows there indeed exists a disturbance in the E region on 13 April. Similar to the abnormal region of electron density, the SNR of orbit no. 030939-1 on 13 April also decreased in the southwestern direction in Fig. 2. This phenomenon may illustrate the decrease in SNR caused by TEC enhancement. Furthermore, this TEC enhancement was probably caused by an earthquake because it shows very intensive conjugate response. However, TEC anomalies caused by geomagnetic storms do not exhibit this kind of phenomenon generally (Zhao et al., 2008).

Figure 6 The spatial distribution of the GPS–TEC map (a, b) and its anomalies (c, d). The GPS–TEC map on 13 April at 06:00 UT (a). The sliding mean of 11 d of background (b). The global anomalies in the GPS–TEC map (c). The regional anomalies around the epicenter of the Yushu earthquake in the GPS–TEC map (d). The purple pentagram indicates the epicenter, and the radius of the black circle is 550 km.

3.2 The possible mechanism of SNR variation revealed by full-wave simulation

In Sect. 3.1, we analyzed the spatial and temporal characteristics of SNR during the five revisit periods before and one revisit period after the Yushu earthquake. It can be found that SNR decreased significantly before the earthquake over the epicenter area of the Yushu earthquake, especially in the southwestern direction. After excluding the influence of geomagnetic storms, we further explored the possible mechanism of SNR abnormal variation in this section. As mentioned in the Sect. 1, the electron density in the lower ionosphere can be disturbed through various mechanisms before earthquakes. The electron density before the Nepal earthquake was obtained from a computer ionosphere tomography method by using GPS data (Kong et al., 2018). Their results shows that the abnormal variation of electron density occurred at a height of 150 km before the Nepal earthquake and the range of variation reached about 30 %. However the electron density hardly changed at a height of 450 km. Marshall et al. (2010) have shown that 60 horizontal discharge pulses of 7 V m⁻¹ near the ground can cause 50 % change in electron density in the lower ionosphere, and 60 horizontal discharge pulses of 10 V m⁻¹ near the ground can even cause a 400 % change in electron density. The variation of electron density in the ionosphere caused by lightning activity and an earthquake can both be explained by one lithosphere–atmosphere–ionosphere coupling mechanism with penetration of the direct current (DC) electric field (Zhou et al., 2017; Kuo et al., 2011). These results provide us a reference for the amplitude of the perturbation of the electron density in the D–E region. Based on these results, the full-wave model was used to simulate the changes in the electric field at satellite altitude excited by ground-based VLF transmitters caused by the enhancement of or decrease in electron density in the lower ionosphere so as to further determine the change law of SNR.

Figure 7 The electron density obtained from COSMIC data on 13 April in the TEC abnormal region.

As mentioned in the introduction, the major VLF wave energy almost lost in the D–E region; after that, the radio waves penetrate the topside of the ionosphere and even magnetosphere with a minor linear reduction because of the mode conversion (Lehtinen and Inan, 2009; Shao et al., 2012). The data of COSMIC also illustrate that the anomaly of electron density not only occurred in the F region (represented by the anomaly of TEC) but also occurred in the D–E region, so the full-wave method (FWM) (Lehtinen and Inan, 2009) was utilized to simulate the electric field between altitudes of 0 and 120 km induced by the NOV transmitter, which is the closest transmitter to the epicenter of the Yushu earthquake. Considering that the study area is much smaller than the radius of the earth, the earth's curvature was neglected in this study. A Cartesian coordinate system was established with x, y in the horizontal plane and z as vertical upward.

Figure 8 The electron density profiles during nighttime. IRI represents the original electron density predicted by the IRI model; IRI+ represents the electron-density-added Gaussian shape perturbation; IRI- represents the electron-density-subtracted Gaussian shape perturbation.

We set a Gaussian shape perturbation at 110 with 20 km standard deviation in the ionosphere. The magnitude of the perturbation was set to a maximum of 1.3 and 4 times both the increase and decrease compared to the original electron density of nighttime (the average electron density above the NOV transmitter on 2–14 April 2010 at 22:00 LT calculated from the IRI-2016 model). The perturbation patterns are shown in Fig. 8 using 4 times the increase and decrease compared to the original electron density as an example. The electron collision frequency is modeled by the exponential-decay law described in Sect. 2.3. The geomagnetic-field intensity and inclination at the location of the NOV transmitter are calculated by the IGRF model.

The electric field only from the ground surface to 120 km has been calculated by the full-wave model because the electromagnetic wave at VLF band will propagate upward as the whistler mode. The group velocities of the upward-radiated whistler mode are almost parallel, and these waves form a narrow-collimated beam, which does not have much lateral spread. The direction of group velocities is determined by the refractive-index surface. The refractive-index surface of the upgoing whistler mode at 120 km is shown in Fig. 9. A ducted propagation is adopted at this L shell (Clilverd et al., 2008), and the VLF wave power is spread in accordance with the divergence of geomagnetic field lines with a linear reduction because of the mode conversion (Lehtinen and Inan, 2009; Shao et al., 2012). The abnormal region of TEC and SNR both that occurred in the southwestern region of the Yushu epicenter could demonstrate that the VLF radio wave propagated in a ducted mode.

Figure 9 The refractive-index surface at 120 km. Red line shows a slice of the refractive-index surface at n_y=0 of the whistler mode, calculated for f=119 kHz at the altitude of h=120 km. Black dash line shows the direction of the geomagnetic field.

The simulated results of the electric field at 120 km height with a different electron density along the magnetic meridian plane within 1000 km area around the transmitter NOV with 11.9 kHz transmitting frequency are shown in Fig. 10. The simulated results are similar when the transmitting frequency is 12.6 and 14.9 kHz. It can be seen that the wave mode interference in the waveguide has been mapped into the ionosphere in the electric field (Lehtinen and Inan, 2009), and the electric field increases when the electron density decreases, and vice versa (Fig. 10a, c). Furthermore, the maximum value of the electric field varying with height is collected to study the influence of the electron disturbance. At nighttime, when i the variation of electron density is smaller, the variation of the electric field is also smaller (Fig. 10b, d). When the electron density increases by 4 times, the maximum electric field decreases about 2 dB at 120 km (see Fig. 10d). The variation is also 2 dB at DEMETER's altitude (660 km) because of the linear reductions (Lehtinen and Inan, 2009; Shao et al., 2012), which implies that the disturbed electric field decreased 20 % compared with the original electric field (Fig. 8b). In the short time interval of a few days before the earthquake, the background noise can be assumed to be stable, so the change in the electric field can reflect the change in SNR. It can be concluded that when the electron density increases by 4 times, the variation of SNR is 20 %. The simulated results illustrate that the variation of electron density in the lower ionosphere before an earthquake is one main factor of causing the abnormal variation of SNR. The more precise SNR variation needs more observation and simulation in the future.

Figure 10 The total electric field excited by the ground-based VLF transmitter NOV with a transmitting frequency of f=11.9 kHz and power of P=500 kW. The total electric field at the altitude of 120 km (a) and the maximum electric field varying with altitude (b) at nighttime when the Gaussian shape disturbance is set to 1.3 times compared with original electron density. The total electric field at the altitude of 120 km (c) and the maximum electric field varying with altitude (d) at nighttime when the Gaussian shape disturbance is set to 4 times compared with original electron density.

4 Discussion

4.1 The possible mechanism on how the earthquake induces the disturbance in the lower ionosphere

Which coupling mechanism is effective to induce electron density anomalies in the D–E layer by earthquakes is still an open question. Molchanov et al. (2006) declared that the lower-ionospheric disturbance is caused by acoustic gravity waves triggered by earthquakes. At present, the coupling mechanism of the electric field proposed by Pulinets (2009) is widely accepted because it has been demonstrated by a series of models (Kuo et al., 2011; Namgaladze et al., 2013; Zhou et al., 2017) and observations (Gousheva et al., 2006, 2008; Li et al., 2017). As for the 2008 Wenchuan Ms 8.0 earthquake in China, Li et al. (2017) reported continuous observations about the anomalous electric field which lasted longer but weaker than the electric field induced by lightning during 1 month before the Wenchuan earthquake, which suggests that the abnormal electric field might be caused by the seismogenic activity of the Wenchuan earthquake. Xu et al. (2011) also found about 2 mV m⁻¹ anomalous electric field in the F2 layer of the ionosphere before the Wenchuan earthquake. Gousheva et al. (2006, 2008) revealed a large number of anomalous electric fields before earthquakes using the INTERCOSMOS satellite. In addition, it is demonstrated that the anomalous electric field induced by an earthquake could change the electron density in the lower ionosphere by Kuo et al. (2011) and Zhou et al. (2017). Such as for the 2015 M 8.1 Nepal earthquake, the electron density variation was well explained by the ground electric-field coupling model established by Zhou et al. (2017).

Figure 11 The Kp and Dst index in April 2010 (a). The SNR distribution on 5 and 6 April with a geomagnetic storm and 13 April (1 d before the Yushu earthquake) (b). Please note that the date format of the heading is year month day, but the date format of the black text is day/month, and the date format of the red and blue text is month day.

4.2 The other factors may induce disturbance in the lower ionosphere

The lightning, geomagnetic storms, and other natural sources may induce disturbance in the lower ionosphere (Marshall et al., 2010; Maurya et al., 2016; Peter et al., 2006; Zigman et al., 2007). As known, the intensive TEC change occurs during geomagnetic storms, and the change in TEC is affected intensively during the main phase of the geomagnetic storm, gradually returning to normal accompanying the recovery phase. To avoid the effect of geomagnetic storms, the data which Kp > 3 and Dst $<-\mathrm{30}$ nT were excluded in this research, and the TEC anomaly detected in Fig. 6 was seen 1 d after the recovery phase of the geomagnetic storm (Fig. 11a). Furthermore, the change pattern of TEC is totally different from the one caused by an earthquake because the TEC anomalies caused by a geomagnetic storm expand from high latitudes to mid latitudes due to thermospheric neutral winds, E × B convection, and so on (Pokhotelov et al., 2008). From Fig. 11b, we can see SNR on the whole orbit are large on 5–7 and 11 April during geomagnetic storms, especially at the higher latitudes. However, the SNR pattern on 13 April is totally different; SNRs on the orbit of 13 April only decrease in the abnormal TEC region. In sum, the TEC anomaly on 13 April should be unconcerned with the geomagnetic storm. A lightning flash is very rare in our research region (only four events from February 2010 to April 2010, which can be seen from the search result of https://lightning.nsstc.nasa.gov/nlisib/nlissearch.pl?coords=?579,18, last access: 27 August 2020), so the effect of lightning could be ignored in this study.

5 Conclusions

In this paper, the SNR of the electric field from a ground-based VLF transmitter observed by the DEMETER satellite was analyzed before and after the 2010 Ms 7.1 Yushu earthquake. The VLF signals from Russian VLF transmitters can be clearly observed at frequencies of 11.9, 12.6, and 14.9 kHz over the epicenter from the electric-field spectrum data. To determine whether the SNR variation is related to the Yushu earthquake, the data in quiet space weather conditions (Kp ≤ 3 and Dst ≥ −30 nT) have been selected during five satellite revisit periods before the earthquake and one revisit period after the earthquake. The result shows that SNR decreased during one revisit period before the Yushu earthquake in all cases. Our analysis of SNR variation also shows that SNR on 13 April is smaller than that on other days over the epicenter; the day-to-day variation of revisit orbit also demonstrates this point, and the decrease in SNR is the most intensive in the southwestern region when we divide the space over the epicenter of the earthquake into four regions. These results are consistent with the TEC anomalies in Fig. 6. In addition, we also analyzed the SNR changes over the epicenter in the same period from 2007 to 2010 as a background map and found that the SNR change trends of one revisit period before the earthquake relative to background time were contrary to those in 2010. The change trend of SNR decreased in 2010 but increased in background time in the first revisit period before the earthquake. The change trend of SNR is the same in other revisit periods both in 2010 and background time. In sum, it can be concluded that the SNR over the epicenter of the Yushu earthquake decreases abnormally in one satellite revisit period before the earthquake, especially in the southwestern region of the earthquake, which is consistent with the observed TEC anomaly before the earthquake. The decrease in SNR before the Yushu earthquake may be due to the enhancement of electron density.

The electron density in the lower ionosphere may change abnormally before an earthquake through some coupling mechanisms. The full-wave simulation result on NOV transmitter, which is the nearest transmitter next to the Yushu earthquake, indicates that the electric field at the altitude of a satellite will change when we add a disturbance of electron density in the lower ionosphere. That is to say that the SNR of the electric field will also change when the background noise is considered to be invariable a few days before the earthquake. The simulated results show SNR does not always decrease before an earthquake like some previous reports show (He et al., 2009; Molchanov et al., 2006; Yao et al., 2013), which depends on the change in electron density. The SNR of the electric field will decrease with the increase in electron density in the lower ionosphere; SNR will increase with the decrease in electron density in the lower ionosphere. It can be concluded that the variation of electron density before earthquakes may be one important factor influencing the variation of SNR.

We will continually explore the law of SNR change and verify the mechanism we proposed with more seismic events by utilizing the newly launched LEO electromagnetic satellite (China Seismo-Electromagnetic Satellite) (Shen et al., 2018; Zhao et al., 2019) in upcoming work.