Selectivity of seismic electric signal (SES) of the 2000 Izu earthquake swarm: a 3D FEM numerical simulation model.

Although seismic electric signal (SES) has been used for short-term prediction of earthquakes, selectivity of SES still remains as one of the mysterious features. As a case study, we made a numerical simulation based on a 3D finite element method (FEM) on the selectivity of SES observed in the case of the 2000 Izu earthquake swarm. Our numerical results indicated that the existence of conductive channel under Niijima island could explain the reported SES selectivity.

Introduction

There are numerous reports on geoelectric changes prior to earthquakes.[1)]–[6)] The so-called seismic electric signal (SES), which is a kind of precursory geoelectric potential changes, has been applied to short-term prediction of earthquakes in Greece.[3)]–[4)],[7)]–[8)] This approach is known as the VAN (Varotsos, Alexopoulos, and Nomicos) method. There have been, however, some active debates on this method.[9)]–[10)] The main concerns include the physical generation mechanism of SES[8)],[11)] and the selectivity phenomenon of SES.[8)],[12)] Selectivity means that SES can be detected only at limited sites, while most of randomly chosen sites are insensitive, and a sensitive site is sensitive only to SES from some specific focal area(s). In this study, we will focus on the SES selectivity.

In order to explain the SES selectivity, some laboratory analog experiments on a geographical scale model and a waveguide model were developed to simulate the propagation of seismic electromagnetic signals.[13)]–[16)] These experimental results showed that the geographical effect such as the distribution of ocean and land may lead to some aspect of the selectivity phenomenon.

Varotsos et al.[17)] proposed a conductive channel model to explain the SES selectivity and presented supporting analytical solutions of Maxwell equations.[18)] The same idea was achieved numerically by Sarlis et al.[19)] using the EM1DSH software, which was developed for calculating either the layered earth model including displacement currents or the 2 layer model with embedded thin sheets but without displacement currents.[20)] Their numerical simulation indicated that there is some amplification effect of electromagnetic signals above the upper end of a conductive channel.[19)] Their work about underground conductive channel did provide an explanation for the SES selectivity. But their numerical simulation can deal with only a 2 layer earth with a thin sheet embedded in the second layer due to the limitation of the EM1DSH software they used.[20)] Thus, whether or not their conclusion holds for a more realistic 3D model deserves further investigation.

The 2000 Izu earthquake swarm started at the end of June, 2000 and lasted for several months.[21)],[22)] The number of events that occurred during June – August, 2000 was more than those of the past 5.5 years in the same region. It had some 7000 M ≥3 shocks and five M ≥6 shocks. The swarm started under Miyakejima island with some volcanic activity on June 26 and rapidly migrated northwestward. The swarm is now considered to have been related with some quasi-horizontal underground magma migration.[23)],[24)] There are reports on the SES activity related to this Izu earthquake swarm,[25)]–[27)] i.e., clear geoelectric disturbances were detected for a period of three months immediately prior to the onset of the swarm. Selectivity of this SES activity was very remarkable because SES was observed only by the Wakago (abbreviated as Wak) measuring dipoles of the Niijima station, while no clear changes were observed by other dipoles of the same and other stations in Izu islands.[25)]–[27)] Strikingly similar observations have been reported in Greece, for example, at the Keratea station close to Athens[28)] and at the Ioannina station in northwestern Greece.[29)]

The aim of this study is to provide a model of the above cited SES selectivity of the Izu earthquake swarm observed at the Niijima station by using a 3D finite element method (FEM) numerical simulation.[30)]

Finite element method (FEM)

The principle of FEM consists of division of the solution domain into small sub-domains, which are called “finite elements”.[30)],[31)] One merit of FEM is its flexibility in simulating a model with arbitrary geometries. Because the reported SES selectivity associated with the Izu earthquake swarm is based on the DC electric potential data, this study will focus only on the case of a DC source, which satisfies the following Poisson’s equation,

where u is electric potential, σ is conductivity, I is DC current source, and δ is Dirac delta function.

Because the total potential u is singular at the source point, which would lead to large error in modeling, a conventional approach is to let u = u 0 + u, where u 0 is the potential due to the source current in a uniform half-space of conductivity σ0, and u is the anomalous potential due to the conductivity anomaly σ, which satisfies σ = σ − σ0.[30)] Then the equivalent variational problem can be described by[30)],[32)]

and

where F(u) is the variational function (energy functional), V is the computation domain, Γ∞ is the external surface of the model boundary, Γ is the boundary including Γ∞ and the air-earth interface, r is the radial vector from the source point, and n is the outward normal direction on the boundary.

With the variational method, instead of solving directly the governing equation, we minimize the above energy functional, F(u), by eq. [3]. There is no source term in eq. [2], thus, the singularity is removed in eq. [2]. The u 0 can be obtained analytically and u can be calculated using FEM.[30)]–[32)]

In this study, we adopted the COMSOL Multi-physics®, a commercial FEM software, to simulate the response to an electric dipole current source in an inhomogeneous medium. As a case study, we made the simulations for Niijima island, at which SES selectivity associated with the 2000 Izu earthquake swarm has been observed.[25)]–[27)]

Validation test

As a validation test of the COMSOL FEM software, we calculated the electric field on the top surface of a three-layer model with a DC dipole source in the 3rd-layer (Fig. 1(a)). The numerical results (circles) obtained from COMSOL are in good consistency with the theoretical ones (black line) (Fig. 1(b)). The maximum relative error is less than 0.5% (Fig. 1(c)).

Fig. 1

Validation test of the COMSOL FEM software. (a) A three-layer model with a DC electric dipole current source. (b) The numerical results (circles) are consistent well with the theoretical ones (black line). (c) The relative error between the numerical and theoretical results.

We also calculated the response of an infinite cylindrical conductive channel lying in a full space medium with a DC dipole source and compared the numerical results with the analytical ones derived by Varotsos et al.[8)],[18)] Again, the numerical results fit well with the analytical ones.

The above tests supported that the COMSOL FEM software can be applied to simulate inhomogeneity models with a DC dipole source.

Numerical model

We assumed a 3D simple layered model as shown in Fig. 2(a). The most important inhomogeneity in this model is the distribution of land and sea. Inhomogeneity within Niijima island may be related with its geological structure. Niijima island is made almost exclusively of rhyolite, except a small outcrop of basalt in its northern tip (Wak area in Fig. 2(b)), the only area where SES was observed.[25)]–[27)] Because of this geological condition, it was suspected naturally if the resistivity difference between rhyolite and basalt may play some important role for the selectivity, and models were examined with varied combinations of their resistivities. More details on the electrical structure of the island being unavailable, our basic 3D model was as simple as shown in Fig. 2. Fig. 2(d) indicates the adopted finite elements (cells), showing cells are much finer for the interested target regions such as Niijima island. After examining the surface inhomogeneity effect, we moved on to a model which contains an underground convective channel as shown in Figs. 2(c) and (d).

Fig. 2

Sketches of the numerical model adopted in our simulation. (a) The 3D model of the calculated region including the ocean-land distribution. The depth of sea water is assumed to be 300 m. The resistivity parameters of sea water and host rock are assumed as 0.5 Ω · m and 2000 Ω · m, respectively. Niijima island is modeled with a 300m thick layer, consisting of rhyolite, with a small basalt exposure in the northern tip (Wakago area). (b) A horizontal view on the ground surface of the 3D model. The two arrows respectively indicate the Wakago (Wak)-Air and the Boe-Air measuring dipoles. SES was observed only in Wakago area where basalt is exposed. (c) An E-W vertical cross-section passing through the hypocenter of the M6.3 earthquake and the vertical conductive channel. The star indicates the location of the current source in most calculations. (d) A 3D view of the cells of the 3D model adopted in FEM calculations.

Surface inhomogeneity effect

In order first to evaluate the possible effect of the inhomogeneity of the surface layer, we assumed an inhomogeneous surface layer with a thickness of 300 m between a homogeneous half-space of air and a homogeneous host rock layer of resistivity 2000 Ω · m (Figs. 2(a)–(c)). A horizontal DC dipole source in E-W direction with a moment of 104 A · m is assumed in the host rock at 10 km depth below the epicenter of July 15 M6.3 earthquake (marked by the star in Figs. 2 and 3) of the Izu earthquake swarm. As stated earlier, the first order inhomogeneity of the surface layer is due to the ocean (resistivity of sea water is taken as 0.5 Ω · m) and Niijima island. Since the island consists of rhyolite and a small outcrop of basalt at Wak (Fig. 2(b)), their electrical property difference would provide the second order inhomogeneity which might be important for the observed selectivity. Different rocks such as rhyolite and basalt have different electrical resistivity, but not in the order of magnitude. The electrical resistivity of the upper crust may vary several orders of magnitude mainly depending on the existence of fluids, such as water, in molecular as well as geological scales.[33)],[34)] Without any more relevant information on Niijima island, we tentatively assumed varied resistivity combinations of basalt and rhyolite as 200/2000, 2000/2000, and 20000/2000. Figure 3 shows the electric field distribution on the ground surface for the above three combinations of resistivity. The results indicates that different resistivity ratios lead to some changes of the distribution of electric field in the northern tip, less (more) resistivity making the Wakago area less (more) sensitive. One can also recognize, in Fig. 3, the effect due to the large contrast of resistivity of sea water and land, which was suggested as one possible explanation of SES selectivity based on some analog experiments.[13)]–[16)] But this ocean-land effect can not be an explanation of the observed selectivity of Niijima island that SES was observed only in the Wakago area at the northern tip of the island.

Fig. 3

Effect due to the surface inhomogeneity of resistivity. The electric field due to a DC dipole source under the star at 10 km depth is given for varied combinations of resistivity of basalt and rhyolite, (a) 200/2000, (b) 2000/2000, and (c) 20000/2000.

Effect of Conductive channel

The effect of underground channel has been proposed as the possible explanation of SES selectivity.[17)]–[19)] We calculated the electric field on the surface for the 3D model with an underground conductive channel placed as shown in Figs. 2(c)–(d). Assumed channel was a vertical cylinder with radius of 0.5 kilometer. As an example, other structures were assumed to be the same as those in the case of Fig. 3(b), namely adopting 2000/2000 for the resistivity contrast of basalt and rhyolite. The results for the vertical conductive channel with resistivity of 20, 200, and 2000 Ω · m are given in Figs. 4(a), (b), and (c), respectively. In order to make the results clearer, we show hereafter the electric field distribution only on the island, since the electric field distribution on the ocean surface is quite similar in all cases. Fig. 4(a) shows that the low resistance (20 Ω · m) of the conductive channel makes the electric field in the northern tip of the island enhanced locally as the observed selectivity requires. The above results suggest that the role of underground conductive channel should be more important than the surface inhomogeneity of resistivity at least for the SES selectivity at Niijima station.

Fig. 4

Effect of underground conductive channel in a model with the same structure as in Fig 3(b). The electric field on the ground surface is given for a vertical conductive channel of a cylinder with radius of 0.5 kilometer below the Wakago site extending from a depth of 300 m to 10 km depth. Varied resistivity is assumed for the conductive channel: (a) 20 Ω · m, (b) 200 Ω · m, and (c) 2000 Ω · m. Lower resistivity of the channel results in enhanced sensitivity of the Wakago site.

We also evaluated the possible effect of the position of the underground conductive channel for the same layer model as Fig. 3(b). Figure 5 shows the results for three example cases in which the resistivity of the channel is assumed 20 Ω · m. If a slanted conductive channel directly connects the source and the bottom of surface layer at 300 m depth below the observing site (e.g., the Wakago site, at which anomalous potential changes are detected[25)]–[27)]), the amplification effect due to the conductive channel is much more significant (Fig. 5(a)). If the channel is buried at deeper position, its effect tends to be weaker as seen in Fig. 5(c).

Fig. 5

Position dependence of the effect of underground conductive channel. The electric field distribution for (a) a slanted conductive channel directly connecting the source to the bottom of the surface layer at 300 m depth below the Wakago site; (b) a vertical channel below the Wakago site extending from the depth of 300 m to 10 km depth; (c) a vertical channel below the Wakago site from a depth of 1 km to 10 km depth. The layer model is the same as in Fig. 3(b) and resistivity of channel is assumed to be 20 Ω · m.

Discussion

As a case study, we have so far simply assumed a DC dipole source at 10 km depth below the epicenter of the nearest July 15 M6.3 earthquake of the Izu earthquake swarm (its epicenter is marked by the star in the figures). In fact, however, the reported SES activity with the same selectivity at the Niijima station had lasted for a few months before the onset of the swarm and while the earthquakes migrated tens of kilometers toward Niijima island before July 15 earthquake possibly along with magma migration.[21)]–[25)] Therefore, it is more likely that the SES activity at Niijima island was related with the whole seismic swarm rather than only with the July 15 earthquake. In order to investigate whether or not the above described results of the simulations hold for other earthquakes of the Izu swarm, we made similar simulations on the same model parameters as those of Fig. 4(a) except setting the DC dipole source at 10 km depth below the epicenter (marked by the star in Fig. 6) of an earthquake that occurred at earlier stage of the swarm (June 26–July 1).[21)],[22)] Considering the role of a conductive channel demonstrated in Fig. 5, we assumed an additional similar but wider horizontal conductive channel (its projection on ground surface is marked by the thick dashed white lines in Fig. 6 with a width of 5 km at 10 km depth between the source and the northern tip of Niijima island). The electric field on the ground surface is given in Fig. 6. One can see that although the Wakago site is at the far end of Niijima island from the source, locally enhanced electric field can be expected at the Wakago site.

Fig. 6

The electric field on the ground surface for a model with the same parameters as in Fig. 4(a) but with an underground conductive channel as described below. The channel with a resistivity of 20 Ω · m is located vertically below the Wakago site from a depth of 300 m to 10 km depth as Fig. 4(a) and then extended horizontally to the current source 10 km below a red star, which corresponds to the epicenter of swarm earthquake at earlier stage of swarm. The horizontal channel has a width of 5 km (the thick dashed white lines indicate its projection on ground surface) to simulate the possible future path of the hypocenter or magma migration. The result indicates the high sensitivity of the Wakago site.

One may argue that the resistivity parameters adopted in the current model simulation are too arbitrary. However, the main concern of this study is to provide some qualitative background of SES selectivity. It has been shown that the first order feature of the electric field on the ground surface is dictated by the ocean-land distribution while the sensitivity test using a varied basalt/rhyolite resistivity contrast showed that less resistivity of basalt in Wakago area makes the area less sensitive to SES as shown in Fig. 3. These effects play minor role compared with the conductive channel in generating the observed SES. In fact, the existence of a conductive channel makes the electric field in the northern tip of Niijima island enhanced locally (Fig. 4).

Our results seem to support that conductive channel plays the key role for SES selectivity, and conductive channel abutting at resistive surface layer near the measuring site produces high amplification (Figs. 4(a) and (b), Figs. 5(a) and (b)), being consistent with the previous numerical results based on 2D model.[19)] Of course, however, the real case should be much more complicated, depending on the details of conductive channel, surface inhomogeneity of resistivity, distance between source and observing site, and so on. It is, therefore, very important to know the detailed underground electrical structure in order to obtain a more realistic explanation of the observed data. Definitely, the current numerical results require further support from electromagnetic sounding data at Niijima island, especially about the possible existence of underground conductive channel. The recent work[35)] on the electrical structure of Kozushima island where SES selectivity, which was much different from that of the Niijima case, was reported[12)] may provide useful information for further understanding of the SES selectivity phenomena. Needless to say that the remaining central problem is the physical origin of the conductive channels.

It should also be mentioned that we discussed only the DC electric dipole source in this study. However, the reported seismic electromagnetic signals vary in a broad band of frequency, from DC to VHF (very-high-frequency). Thus, the investigation of the frequency effect of the source would be an interesting topic, which will be discussed in a separated work.