The numerical simulation on ionospheric perturbations in electric field before large earthquakes

Many observational results have shown electromagnetic abnormality in the ionosphere before large earthquakes. The theoretical simulation can help us to understand the internal mechanism of these anomalous electromagnetic signals resulted from seismic regions. In this paper, the horizontal and vertical components of electric and magnetic field at the topside ionosphere are simulated by using the full wave method that is based on an improved transfer matrix method in the lossy anisotropic horizontally stratified ionosphere. Taken account into two earthquakes with electric field perturbations recorded by the DEMETER satellite, the numerical results reveal that the propagation and penetration of ULF (ultra-low-frequency) electromagnetic waves into the ionosphere is related to the spatial distribution of electron and ion densities at different time and locations, in which the ion density has less effect than electron density on the field intensity. Compared with different frequency signals, the minimum values of electric and magnetic field excited by earthquakes can be detected by satellite in current detection capability have also been calculated, and the lower frequency wave can be detected easier.


Introduction
ULF (ultra-low-frequency) electromagnetic anomalies (0-15 Hz) have been widely detected by ground-based stations and satellite observations before many earthquakes.The magnetic field at 0.01-5 Hz increased before the Ms 7.1 Loma Prieta earthquake (Fraser-Smith et al., 1990).Electric field variations of 100-300 mV km −1 began in 3-16 days prior to some earthquakes (Myachkin et al., 1972).An anomalous increase of 3-7 mV m −1 in the vertical component of the quasi-static electric field was observed on board INTERCOSMOS-BULGARIA-1300 satellite in the equatorial ionosphere over the seismic region about 15 min before a M 4.8 earthquake (Chmyrev et al., 1989).Gousheva et al. (2008) found abnormal increase up to 2-15 mV m −1 in the horizontal and vertical components in the quasi-static electric field in the near-equatorial, low, middle, and polar latitude ionosphere before seismic activities.Zeren et al. (2012) also found the electromagnetic perturbations 4 days before Yushu earthquake.Zhang et al. (2014) studied the ULF electric field (DC-15 Hz) observed by the DEMETER satellite (h = 660 km), and found that the components of ULF electric field perturbations were with opposite variation shape in E x and E z , as well as that the increase of electron density and O + density were detected with the disturbances in electric field together.
As for the seismo-ionospheric coupling theory, Nagano et al. (1975Nagano et al. ( , 1993) ) had constructed a model using the full wave method to study the VLF electromagnetic waves (3-30 kHz) penetrating into the ionosphere.The waves excited by earthquakes are mainly at a ULF band, so the effect of the ion cannot be ignored as that at a VLF band.In this paper, a model considering the effect of ions is constructed based on the VLF wave propagation code of Nagano et al. (1975) to study the effects with different backgrounds and variations of the ionospheric ion and electron densities on the electromagnetic field at the topside ionosphere when a ULF wave penetrates from the bottom of the ionosphere.And two earthquakes in equatorial and mid latitudes studied by Zhang et al. (2014) are selected to simulate the abnormal characteristic of electromagnetic field change due to different background and variation of plasma parameters, and the results are compared with the actual observations recorded by the DEME-TER satellite.There have been opposite views about if the electric fields excited by earthquake can penetrate into the ionosphere.Some has concluded that electromagnetic waves can penetrate into the ionosphere (Kim et al., 1994;Pulinets et al., 2000;Sorokin et al., 2001Sorokin et al., , 2007)), but the others hold opposite views, that the electromagnetic energy released by earthquakes is too weak to be observed at satellite altitudes (Densisenko et al., 2008;Ampferer et al., 2010).Therefore, based on the current detection ability of satellites, the minimum electric field at the bottom of the ionosphere is also calculated by our new model in this study.
In Sect. 2 the method is introduced to calculate the electromagnetic field of ULF-wave propagating in the horizontally stratified anisotropic ionosphere.In Sect.3, simulations are done to test and verify the ULF model, and then it is applied in real earthquake research.Finally, the conclusions drawn from the numerical simulation are summarized in Sect. 4.

Method
The electromagnetic field can be got by solving the following Maxwell equations in the anisotropic ionosphere, in which ε 0 , µ 0 , I denote the dielectric constant, magnetic permeability of free space and the unit matrix, respectively.
To solve these equations, the coordinate system is set as followed.The Y axis points to the geomagnetic north direction, and the X axis points to the east.The Z axis is perpendicular to the X-Y plane.The Cartesian coordinate system is shown as Fig. 1.The plane wave is in the X-Z plane with incident angle θ i .The cosines direction of oblique geomagnetic field H 0 is (l, m, n).The ionosphere is idealized as a horizontally stratified anisotropic media in Z direction.
The dielectric susceptibility matrix M in the ionosphere is expressed like Eq. ( 2) by Yeh and Liu (1972).
The subscript α represents the different types of particles in the ionosphere including electron, O + , He + , H + and so on.
, where ν α is the collision frequency of a specific type of particles.
N α e 2 α m α ε 0 represents plasma angular frequency of a specific particle, and m α , e α and N α are its mass, electron charge and density.
represents the plasma gyro angular frequency of a specific particle.
Considering that the ionosphere parameter varies quickly within the scope of one wavelength, a horizontally stratified ionosphere is assumed and the space-time factor of each field component of the incident wave is given as where k 0 is the propagation constant in free space, S = sin θ i , C = cos θ i .
After eliminating the E z and H z components from Eq. (1) using Snell's law, an elegant form of Maxwell's equation is obtained (Budden, 1961): where e is the state vector of each layer composed by horizontal components of electric and magnetic fields T is the state matrix of 4 × 4 in each layer, while the S, C and M in matrix T are expressed as above.
The four eigenvalues and eigenvectors of matrix T represent the index and amplitude of four up-going and down-going waves in the ionosphere, respectively, as shown in Fig. 2. Transfer matrix method is used from the bottom layer to the top to get the field in each layer in Eq. ( 5).where K s denotes the propagation matrix in each layer; A and B are unknown coefficients.These coefficients will be determined by applying the incident field conditions at z = z 1 .
To avoid the numerical swamping resulted from the loss of medium, the orthogonal scale down process developed by Nagano et al. (1975) is also used in this study.Because Nagano's code is only used to solve the electromagnetic field of the VLF wave, the matrix M is simplified as expression (6).In this matrix, only the electron is considered.Compared with the VLF wave frequency, the gyro-frequency of the ion is very small, so it can be ignored (Lu et al., 2008).But at a ULF band, where gyrofrequency of the ion is closed to the wave frequency, the effect of the ions are so vital that they cannot be neglected.The expression ( 6) is not reasonable in a ULF band, therefore we modify the matrix M of Nagano's model from expression (6) to (2) with ions taken into consideration.
3 Simulation of electromagnetic waves penetrating into ionosphere

Calculation parameters
Below 50 km the model has been set to free space, whereas above 50 km the model is set to horizontally stratified anisotropic ionosphere.The intensity of incident wave from free space is set to one unit, by which all the results are normalized.To compute the electromagnetic field at the topside ionosphere excited by the ULF wave by equations described in Sect.2, we need to get the basic ionospheric parameters such as ion densities, electron density, total collision frequency, as well as the magnitude and direction of geomagnetic field.Generally, the International Reference Ionosphere model (Bilitza, 2001) has been used to compute the N e (electron density) and N i (ion density) profile from 50 to 700 km.The model of collision frequency between 80 and 300 km uses the model of Cummer (2000), and Helliwel (1965) Model has been applied from 300 to 700 km.The geomagnetic field is computed using the International Geomagnetic Reference Field (Barton, 1997).It has been observed that the ion density (N i ) and electron density (N e ) rose before earthquakes (Zhang et al., 2014).In this study, we choose two cases on 12 November 2008 and 28 February 2010 when the increase of N e and N i was recorded (Zhang et al., 2014) before two earthquakes on 16 November 2008 and 5 March 2010, and compute N i and N e profiles by multiplying increase factors on the profiles given by the IRI (International Reference Ionosphere) model.Figure 3a, b illustrates the N e and N i profiles at nighttime during those 2 days.The geomagnetic field strength and inclination are computed according to the occurrence time and location of anomalies mentioned by Zhang et al. (2014), which is given in Table 1.
Other parameters like the incident angle of the radio wave is set at 45 • , and the thickness of every layer is set at 3 km in the numerical simulations of this paper.

Verification of ULF propagation model
Our ULF model, which is constructed to calculate the electromagnetic field distribution excited by the ULF wave penetrating from the bottom of the ionosphere, is different than that at a VLF band, so it needs to be verified.
First two simulations are done to test whether ULF model works.The four horizontal components of electric and magnetic field are calculated at two different frequencies, 10 kHz (VLF) and 10 Hz (ULF).At 10 kHz, the frequency is so high that the ions can be ignored, which means the calculated results must be the same, whether using the VLF model or the ULF model.Figure 4 shows the simulated results at 10 kHz by using Nagano's VLF model (Nagano et al., 1975), which are consistent with the results from our new ULF model.It can be seen that the solid lines (ULF model) totally super-  Another simulation is done to test the correctness of our new ULF model.The penetrating processes of 10 and 100 Hz electromagnetic waves are simulated separately using the same ionospheric profiles at (47.0 • S, 91 • W) on 28 February 2010 shown in Fig. 3a.For 10 Hz waves, both "penetrating" (right-handed polarization) and "non-penetrating" (left-handed polarization) modes are able to propagate into the topside ionosphere as shown in Fig. 6a and the wave energy is mainly stored in the magnetic field.However the nonpenetrating mode for 100 Hz wave is strongly damped in the ionosphere (Fig. 6b).It can be concluded that the higher frequency wave experiences severer absorption in the D region of the ionosphere than the lower frequency wave as expected (Helliwell and Bleier, 1965), which coincides with the result of Bortnik and Bleier (2004).

The electric field variations before earthquakes
To illustrate the relationship among electric field and plasma parameters, the intensities of E x and E z at 10 Hz are calculated with four different density profiles, including normal  N e with N i , normal N e with 10 % higher N i , normal N e with 100 % higher N i , and 5 % higher N e with normal N i as listed in Table 2.The normal profiles are calculated by IRI model on 28 February 2010 as described in Sect.3.1, and the other parameters needed in numerical simulation are also given in Sect.3.1.
In Table 2, the simulated results in E x and E z of normal N e with 10 % higher N i is close to that with normal N e and N i , but the attenuation of electric field become larger when normal N e with 100 % higher N i and normal N i with 5 % higher N e are used, which means the effect of ion density on the propagation of the ULF wave is smaller than that of electron density.
As mentioned in Sect.3.1, Zhang et al. (2014) has observed N i and N e increased before two earthquakes located in the middle latitude and equatorial regions separately dur- ing nighttime.And the abnormal time on 28 February 2010 is 7 days before the earthquake taking place on 5 March 2010 at the middle latitude region (47.0 • S, 91 • W), and 12 November 2008 is 4 days before the earthquake on 16 November 2008 located at the equatorial region (0.2 • S, 119 • E).In this simulation the N e and N i profiles around the time and location of the two earthquakes have been calculated with IRI model as normal densities.And an increase factor of 5 % higher N e and 10 % higher N i are both added to normal IRI profiles to describe the observation results that the increase of electron density and ion density before the two earthquakes (Zhang et al., 2014).
The E x and E z components of electric field at 10 Hz and 1 Hz are computed, respectively.shows the results at 1 Hz.It can be seen that the electric field at 1 Hz is larger than at 10 Hz with normal N i and N e , which means that when the frequency is lower; the propagation loss is smaller for ULF wave in the ionosphere (Tables 3 and 4).Another characteristic is that the attenuation of electromagnetic field is smaller in the middle latitude region (20100305) when compared the result in equatorial region (20081116 from Table 3).That means the electromagnetic response of ULF radio wave is much easier to be detected in the middle latitude region.Compared the results by normal and higher N e , the attenuation will strengthen at middle latitude (20100305 in Tables 3 and 4), either with 10 Hz or 1 Hz waves when the electron density increases.However at the equatorial area, the attenuation becomes weak with higher N e (20081116 in Tables 3) by 10 Hz wave.
From the above results in Tables 3 and 4, the largest attenuation of the intensity of electric field is almost 60 dB from the bottom of the ionosphere to the altitude of DEMETER satellite (20081116 in Table 4 with normal higher N e and N i ).The sensitivities of the electric field antenna onboard DEME-TER satellite are 0.1 µV m −1 Hz −1/2 at 100 Hz (Berthelier et al., 2006) which means when the electric field of incident wave exceeds 1 V m −1 at the altitude of 50 km, the disturbance can be detected by DEMETER satellite.Xiong (2001) has pointed out that when the seismic charge density is 10 −4 -10 −5 C m −2 , the electric field can reach 10 3 -10 4 V m −1 at 50 km altitude in the late pre-seismic stages in a 100 km 2 region.It seems that the satellite can detect the electromagnetic field of a ULF band excited by seismic activity with the current ability of detection.

Discussion and conclusion
It is known that the gyrofrequency of an ion is comparable to the wave frequency at a ULF band, so the effect of an ion density cannot be ignored in ULF wave propagation model.In this paper, transfer matrix method considering the effect of ion density has been used to construct a new ULF wave model in the ionosphere to study the effects of plasma parameters (N e and N i ) on the propagating waves at different frequency in a ULF band.N e and N i profiles at different latitudes have been chosen to simulate the variation of electric field related to two earthquakes.In this study, the propagation of waves from the ground to the bottom of the ionosphere has not been included, and it needs to be improved in future.Regardless, some conclusions can be drawn as follows.
Firstly, since the variations of ion density actually influence the electric field at satellite altitude in a ULF band, the past VLF model is not applicable in this situation, therefore a new ULF model has been developed in this research by considering the ion parameters.The calculation results in the same VLF band by past and present models, which validates the correctness of this new one.
Secondly, the effect of electron density on the electromagnetic field is larger than ion density, which means electron density measurement with higher accuracy is demanded.
Thirdly, the attenuation of ULF wave is smaller with lower frequency which means the electromagnetic response of lower frequency wave can be detected easier in the topside ionosphere.
Fourthly, the increase of electron density in the ionosphere always make the strength of electric field reduce at the topside ionosphere, either at 10 or 1 Hz.The attenuation of electromagnetic field is smaller over middle latitude region when compared the result in equatorial region.It illustrates that the ULF electromagnetic wave can penetrate into the topside ionosphere much easier over middle latitude region with higher geomagnetic dip.
Last but not least, the minimum amplitude of incident electric field which can be directly observed by DEMETER satellite has been estimated.Although it is just a rough result, it still can be deduced that to detect the electromagnetic response with satellite in the ionosphere is a effective tool for finding earthquake precursors.

Figure 1 .
Figure 1.The Cartesian coordinate system (k indicates the incident electromagnetic wave at X-Z plane; H 0 is the Earth's magnetic field).

Figure 2 .
Figure 2. The propagating electromagnetic waves in the horizontally stratified media (R means the right-hand polarized wave and L represents the left-hand polarized wave)

Figure 3 .
Figure 3. Ion and Electron densities profiles on 28 February 2010 (a) and 12 November 2008 (b) at local time 22:00 before two earthquakes (their geographical latitude and longitude are written on the top of each image).

Figure 4 .
Figure 4. Intensity of four horizontal components of electric and magnetic field varying with altitude at 10 kHz calculated with Nagano's VLF model (dashed line) and our ULF model (solid line) (blue lines represent E x , green ones for E y , red being H x , and purple for H y ).

Figure 5 .
Figure 5. Intensity of four horizontal components of electric and magnetic field varying with altitude at 10 kHz (solid line) and 10 Hz line) calculated with our ULF model(blue lines represent E x ; green lines represent E y ; red lines represent H x ; purple lines represent H y ).

Figure 6 .
Figure 6.Intensity of four horizontal components of electric and magnetic field varying with altitude.Solid lines show intensity of the non-penetrating mode wave, and dashed lines show intensity of the penetrating mode wave at (a) 10 Hz and (b) 100 Hz.Blue lines represent E x , green lines represent E y , red lines represent H x , purple lines represent H y .

Table 2 .
Simulation results of E x and E z (660 km altitude) with different N e and N i profile at 10 Hz. + N i N e + 10 % higher N i N e + 100 % higher N i 5 % higher N e + N i

Table 3 .
Simulation result of E x and E z (660 km) at 10 Hz during earthquake and before earthquakes.

Table 4 .
Simulation result of E x and E z (660 km) at 1 Hz during earthquake and before earthquakes.
i + N e 5 % higher N e + N i N e + N i 5 % higher N e + N i