Volume 123, Issue 17 p. 8901-8914
Research Article
Free Access

Earthquake Lights: Mechanism of Electrical Coupling of Earth's Crust to the Lower Atmosphere

Jaroslav Jánský

Corresponding Author

Jaroslav Jánský

Communications and Space Sciences Laboratory, Department of Electrical Engineering, School of Electrical Engineering and Computer Science, The Pennsylvania State University, University Park, PA, USA

Correspondence to: J. Jánský and V. P. Pasko,

[email protected];

[email protected];

[email protected]

Search for more papers by this author
Victor P. Pasko

Corresponding Author

Victor P. Pasko

Communications and Space Sciences Laboratory, Department of Electrical Engineering, School of Electrical Engineering and Computer Science, The Pennsylvania State University, University Park, PA, USA

Correspondence to: J. Jánský and V. P. Pasko,

[email protected];

[email protected];

[email protected]

Search for more papers by this author
First published: 28 June 2018
Citations: 6

Abstract

Earthquake lights are an atmospheric luminous phenomenon occurring before, during, and occasionally after earthquakes. In the present work we assume that strong electric currents are produced inside the Earth during earthquakes and model the current and field distribution inside and above the Earth to find the conditions where electrical discharges can appear. The simulation results quantify the effects of various configurations of current dipoles inside the Earth. The results are supported by approximate formulations allowing effective solution of the same problems using an analytical theory. It is shown that a large-scale dipole with poles located in Earth's crust at 5 and 15 km beneath the Earth's surface requires energy significantly higher than the total seismic wave energy in major earthquakes. The more likely setup to produce earthquake lights is found to be when the upper pole of the source current dipole is shifted close to the Earth's surface, in particular, when locations of current dipoles are tens of meters beneath Earth's surface.

Key Points

  • Earthquake lights are studied based on enhancement of electric field near Earth's surface leading to conventional lightning-like gas discharges
  • The most likely setup is found to be when the upper pole of the source current dipole is located close to the Earth's surface
  • Large-scale dipole located 5 and 15 km beneath Earth's surface requires energy significantly exceeding that available in major earthquakes

1 Introduction

Earthquake lights (EQLs) are an atmospheric luminous phenomenon occurring before, during, and occasionally after earthquakes and lasting from a fraction of a second to a few minutes (Derr, 1973; Heraud & Lira, 2011; St-Laurent et al., 2006). Theriault et al. (2014) performed a study to determine the parameters of the best documented earthquakes that were preceded, accompanied, or followed by luminosities. This study provides the best up-to-date overview of well-documented EQLs and corresponding earthquake characteristics. It was shown that EQLs can be split based on the time of appearance in two major groups: (i) pre-earthquake EQLs, which are observed from a few seconds to up to 3–4 weeks prior to seismic events, and (ii) coseismic EQLs, which are observed during the rupture event and the passage of the seismic waves. The coseismic EQLs are most frequently observed as either short flashes of light shooting high up into the air, diffuse, but relatively bright, atmospheric glows, or flame-like luminosities coming out of the ground (Theriault et al., 2014). The above mentioned descriptions of EQLs resemble upward lighting flashes (e.g., Rakov & Uman, 2003, Chapter 6) and corona discharges on ground (e.g., Rakov & Uman, 2003, p. 78) observed during thunderstorms, which are well-documented examples of the luminous phenomenon and can help in understanding EQLs. The discharges beneath thunderstorms appear in the areas with increased electric field; for example, corona around trees and other sharp objects appears at the ambient fields higher than 5 × 103 V/m (Standler & Winn, 1979).

To support the existence of high fields above the ground during earthquakes, the measurements of corona currents on an 8-m tall probe were reported by Kamogawa et al. (2004). The authors observed an increase of corona current during earthquake for several minutes (Kamogawa et al., 2004, Figure 4), which indicated the increase of the ambient atmospheric electric field. The current was compared with current under the thunderstorm conditions, and a similar behavior was observed suggesting that fair weather electric field ∼100 V/m (Williams, 2009) was increased by at least 2 orders of magnitude, analogically to thunderstorms. Additional observations of quasi-static electric field of similar amplitude above the ground are reported in Sorokin and Ruzhin (2015, section 2, and references therein). The electric field perturbations related to earthquakes are also measured inside the Earth. The related results are documented in the seismoelectromagnetics review papers (Hayakawa & Hobara, 2010; Uyeda et al., 2009). The horizontal electric field perturbation inside the Earth is caused by telluric currents during earthquakes (Uyeda et al., 2009). The measurements are usually performed at distances of 10–100 km from an epicenter, and measured fields consist of one pulse with duration of 1 s (1 s is also the time resolution of measurements; Varotsos et al., 2007, Figure1); peak values range from 10−6 to 10−4 V/m (Orihara et al., 2012; Varotsos et al., 2007). These electric signals were observed before the earthquake. It is important to note that coseismic electric signals are also present, and they are observed at the time of arrival of the seismic waves and not at the time of earthquake origin (Uyeda et al., 2009, section 4.1, and references therein).

Details are in the caption following the image
Schematics of the computational domain. Panel (a) represents schematically (drawn not to scale) the full view of the Earth, the top ionospheric boundary, and the zero potential remote boundary. Panel (b) represents zoom-in view of the volume close to the axis θ = 0° where the source current dipole Is is introduced with poles at altitudes h+ and h. For such dipole the positive charge is induced on the Earth's surface above the dipole, while the negative charge is induced on the Earth's surface everywhere else. The peak electric field E0 is found directly above the dipole on the Earth's surface. The quantity depicted as σatm corresponds to the conductivity of the atmosphere, and σEarth corresponds to the conductivity of the Earth.

Several modeling studies have been performed to explain the seismoelectromagnetics' observations. Molchanov et al. (1995) assumed a current source inside the seismic source, 10 km beneath the Earth's surface. The authors obtained agreement with the ultralow-frequency (ULF) magnetic and electric signals observed on the ground. Bortnik et al. (2010) assumed the current source also 10 km beneath the Earth's surface. The authors showed that to reproduce their reference earthquake observation, the expected telluric currents fall in the range ∼104–105 A. Both of the above mentioned studies use the current dipole as their current source. Kuo et al. (2014) schematically show in their Figure 3 a dipole representing source currents in the lithosphere and currents flowing in the atmosphere. The origin of source currents is not known at the moment, but several theories are presented in the seismoelectromagnetics literature (Hayakawa & Hobara, 2010; Uyeda et al., 2009). In particular, Freund et al. (2006) developed a laboratory experiment to demonstrate that transient electric currents can flow out of a mechanically stressed rock. The current flowing out of the stressed rock volume is carried by p-holes, that is, defect electrons on the oxygen anion sublattice associated with energy states at the top of the valence band. Because of their wide spectrum of lifetimes, the number of positive holes activated per unit volume of rock and, hence, the outflow currents increases steeply with the rate at which stresses are applied, reaching values as high as 109Akm−3 (Scoville et al., 2015). The coactivated electrons are trapped in shallow energy levels below the edge of the valence band and cannot flow out of the stressed rocks. Freund et al. (2009) have shown that charge buildup at the rock-air interfaces can lead to an enhancement of electric field around sharp points and therefore to air ionization. Their experiments reported current density of 0.1 nA/cm2. The buildup of stress within the Earth's crust prior to major earthquakes may lead to similar manifestations at the Earth's surface.

The motivation of the present work is to find out the electrical conditions for which the EQLs, in the form of electrical discharges, appear. First, we assume the current sources located ∼10 km beneath the Earth's surface to see if the quasi-static electric field can reach values allowing for a discharge ignition above the ground. We note that electromagnetic waves produced by this source can be simulated and detected (Bortnik et al., 2010), but the quasi-static electric field is necessary to ignite the electrical discharge. Second, we assume small-scale current sources ∼100 m beneath the Earth's surface to see whether the local source, possibly related to seismic wave propagation, can be responsible for EQLs. Additionally, we also compare required electrical energy of current source with total seismic wave energy. Finally, we present a comparison of the horizontal electric field just beneath the Earth's surface at a distance of 100 km from the source with experimental measurements of telluric currents (Varotsos et al., 2007).

To achieve the above mentioned goals, we employ the model developed for description of global electric circuit (GEC) (Jánský & Pasko, 2014). This model uses a current source as an input and provides solution of the electric field in the atmosphere and in the Earth. We vary the current source to achieve the field above the ground necessary for discharge ignition. Then, the GEC model is directly used for evaluation of the electrical energy from the first principles. As there is no available information on the electrical energy in earthquakes, our results are then compared with the total seismic wave energy We in units of Joules (J) released in an earthquake logWe=5.24 + 1.44M (Båth, 1966, equation (11.9)), where M corresponds to the magnitude of earthquake on Richter scale. We assume that for the feasible physical scenarios the electrical energy should be smaller than the total seismic wave energy.

The results of present work (section 3), the source currents, and the required electrical energy can serve as references for determining feasibility of different origins of the source current during earthquakes. The simulation results are accompanied with an approximate analytical formulae (section 7) providing means to make quantitative estimates, including those for cases considered in the present work.

2 Model Formulation

We use a time-dependent axially symmetric spherical model of the electrical behavior of the GEC described in Jánský and Pasko (2014). Briefly, the model is based on the continuity equation for charge density ρ
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0001(1)
and the Poisson's equation for the electric potential ϕ
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0002(2)
where urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0003 is the electric field, σ is the conductivity of the media, Scur is the current source term, and ϵ0 is the permittivity of free space.

A schematic of the simulation domain is shown in Figure 1a. The model employs spherical coordinate system (r,θ,φ) with the origin positioned at the center of the Earth at r = 0 km. The model domain includes the highly conducting Earth, a relatively thin layer representing weakly conducting atmosphere and the ionosphere, and a large free space domain bounded by a remote spherical shell representing a physical condition of zero potential ϕ = 0 at infinity. The radius of the remote spherical shell is rout=937rEarth=6.0 × 106 km. Equations 1 and 2 are discretized using a finite volume scheme on a structured mesh, with the mesh refined below the Earth's surface to 10 m in both radial direction and horizontal directions.

A typical range of conductivity of Earth's materials is 10−5–102 S/m (Palacky, 1988, Figure 2). The Earth's crust vertical conductivity profiles are documented in Schwarz, (1990, Figures 1 and 7) and exhibit the variability from 10−5 to 100 S/m. The conductivity profiles are also dependent on geographic longitude and latitude (Palacky, 1988). In the present work the conductivity inside the Earth is set to urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0004 to represent the minimum conductivity of the Earth's crust as it provides the best case scenario for existence of EQLs. The discussion on influence of conductivity is provided in section 8 using an approximate analytical formulation.

Details are in the caption following the image
Two-dimensional view of (a) magnitude of the electric field and (c) charge density around the source currents. The volume above the Earth's surface with the electric field higher than 105 V/m is considered as the place where earthquake lights can appear. The source currents are placed at 5 and 15 km beneath the Earth's surface. The time evolution of the peak electric field E0 above the Earth's surface is shown in panel (b). The time evolution of the external work of the earthquake is shown in panel (d) and compared with the total seismic wave energy of earthquake (M = 7.5) 1016 J. Panel (b) includes the field evolution for the case with relative permittivity of the Earth ϵr=60. Panels (b) and (d) report the field and work time evolutions for pulsed source currents with the sine function over half period τ = 1 s. Panel (d) contains the energy required for the horizontally large source currents based on the work of Kuo et al. (2014).
In the air above the Earth surface an exponential approximation of the conductivity is used (Dejnakarintra & Park, 1974):
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0005(3)
where h is the altitude and l = 6 km is the altitude scaling distance. Then rTIB=rEarth+hTIB, where hTIB=146 km, describes the radial distance from the Earth surface to the top of the model ionosphere. Above rTIB, in the area considered as free space, the conductivity is set to 0. The zero conductivity of free space surrounding GEC keeps charge confined to GEC. The exponential conductivity profile is considered as a very good approximation below 70 km, where ion conductivity prevails (Dejnakarintra & Park, 1974, Figure 4). The high altitude of the top ionospheric boundary hTIB is used to capture principal dynamics of the real system by employing exponential conductivity profile and simple boundary conditions (Jánský & Pasko, 2014).

The relative permittivity of rocks inside the Earth can exhibit permittivity between 5 and 60 for the low frequencies ∼100 Hz, and it depends strongly on the water content (Fuller & Wait, 1976). We do not study this large range of permittivity as the results relevant to the EQLs are independent of the value of permittivity. The reason for this is as follows. The steady state inside the conducting Earth is achieved on a dielectric relaxation timescale ϵ/σec, which is ∼10−6 s for ϵr=1 and ∼10−4 s for ϵr=100. Both dielectric relaxation times are shorter than time of EQLs. Once the Earth reaches the steady state, the results of the present model are a function of the conductivity and independent of the permittivity of Earth. It can be seen in the approximate formulations presented in sections 3 and 7 that they are all independent of the relative permittivity of Earth. In the present work the results are shown for Earth with relative permittivity ϵr=1. An additional test simulation was run, and the independence of results relevant to EQLs for ϵr=60 was confirmed and is presented in Figure 2b.

The zoom-in view of the volume close to the axis θ = 0°, where the source current dipole Is is introduced at altitudes h+ and h, is shown in Figure 1b. The values of h+ and h will be mentioned and justified in section 3 for each earthquake scenario. For a dipole with the upper positive pole, the positive charge is induced on Earth's surface above the dipole, while the negative charge is induced on the Earth's surface everywhere else. One of the parameters reported in our study is the peak electric field E0 directly above the dipole on the Earth's surface. Its location is shown also in Figure 1b.

The quantitative information about the region of charge separation inside the Earth is not available. The approximation of this region as a current source dipole is the basic assumption intended to employ the simplest physical formulation. We note that any complex charge separation inside the Earth can be decomposed into a sequence of dipoles analogically to GEC models where any thundercloud current could be represented using monopoles (Jánský & Pasko, 2014). For example, if the charge separation/current generator is moving along the Earth's surface with seismic waves, the current generator could be described as a sequence of current dipoles moving along the Earth's surface with the same speed.

For the electrical discharge to appear in the air at the ground pressure the field has to exceed the critical field for air breakdown ≃3 × 106V/m (e.g., Morrow & Lowke, 1997). Around local sharp points (i.e., vegetation) the field is enhanced and a lower ambient field is needed for ignition of the electric discharge. For example, the corona around trees and other sharp objects beneath thunderstorms appears at fields higher than 5 × 103 V/m (Standler & Winn, 1979). The EQLs are expected to appear at intermediate fields between the above mentioned values. For quantitative calculations presented below we use a representative reference field 105 V/m.

To evaluate if the considered source current is realistic, we compare the total seismic wave energy in an earthquake with the external energy input required in our simulation. We assume that the electrical energy input from earthquake should be much smaller than the total seismic wave energy. As already noted in section 1 the seismic wave energy in units of Joules (J) is estimated as logWe=5.24 + 1.44M (Båth, 1966, equation (11.9)). For comparison the reference energy 1016 J, corresponding to magnitude M = 7.5 on Richter scale, is used. The external energy input from an earthquake is evaluated as work Ws, which earthquake needs to perform in order to support generation of the source currents. The power supplied by earthquake in volume V of the whole computational domain is urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0006, where the current density urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0007 represents external electrical sources (Balanis, 2012, equation (1-59a)). In the identity based on the divergence theorem urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0008, the left-hand side is 0 because ϕ = 0 on the boundary of the computational domain. The source current density urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0009 is linked with the source current term urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0010 (Browning et al., 1987, equation (2)). Then the external energy input of an earthquake from the start at t = 0 to time t can be written as
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0011(4)
where the integration is performed over the volume Vol with the nonzero current source term Scur.

3 Numerical Model Results

This section presents the results for the electric field, the charge density, and the dissipated energy for several cases of source current setup. First, the scenario representing the current source inside the seismic source is approximated by the current dipole at h=−15 km and h+=−5 km. Second, we assume that the current is externally driven from the seismic source to the Earth's surface represented by the current dipole at h=−15 km and h+=0 km. Third, we examine a small-scale scenario with h=−150 m and h+=−50 m that can be associated with the seismic waves instead of a seismic source.

3.1 Case With h=−15 km and h+=−5 km

Theriault et al. (2014) present a simplified tectonic model to explain the appearance of the EQLs. In their model the charges flow from the stressed volume in the Earth's crust (Theriault et al., 2014, Figures 8 and 9). In this section the representative values for the source current locations are therefore selected as h=−15 km and h+=−5 km corresponding, respectively, to the lower crust and the upper crust. The source currents are homogeneously distributed in the cylinders centered at h and h+ with the radius of 1 km and the height of approximately 3 km.

To obtain the field above the ground E0 higher than the EQLs reference threshold 105 V/m, the source current has to be set to Is=109A. A two-dimensional view of the magnitude of the electric field is shown in Figure 2a. The electric field is taken at time moment t = 1 s representing the typical timescale of EQLs. The volume where the electric field is higher than 105 V/m is enclosed by the black line, and it has the radius of about 10 km and the height of 5 km above the Earth's surface. The peak field E0 in the atmosphere is located right above the dipole at the Earth's surface. In order to illustrate the field geometry, the arrows are added to represent the electric field direction. It is seen that in general the field is directed from the positive pole to the negative pole as expected for a dipole in free space. The difference is that the high conductivity of the Earth maintains the electric current inside the Earth and therefore the field lines that would direct the current outside of the Earth are constrained inside the Earth just below the surface. Only a small fraction of the current (2.2A) escapes to the atmosphere. Note that the magnitudes of the electric field below and above the surface shown in Figure 2a are comparable. Therefore, the ratio of currents below and above the surface is similar to the ratio of conductivities 5 × 10−14/10−5=5 × 10−9.

The field inside the Earth just beneath the surface at 100 km distance from the dipole is 50 V/m. This field is much higher than the measured values 10−6–10−4 V/m (Orihara et al., 2012; Varotsos et al., 2007). If the current source would be selected as 104A as used in (Bortnik et al., 2010), the field would decrease proportionally to the current to a value 5 × 10−4 V/m and would be in good agreement with the experimental measurements. This may suggest that the source current at the seismic origin, which may be responsible for the telluric currents, is not responsible for the discharges above the ground. However, Varotsos et al. (1998) demonstrated that the field profiles created from the dipole in the Earth can change dramatically based on the assumed conductivity distribution. Therefore, the comparison of the computed fields with the measured fields provides an important reference but does not cover all the possible parametric space.

The time evolution of the peak electric field in the atmosphere E0 is shown in Figure 2b. The initial conditions correspond to a zero charge in the domain and therefore to the zero field. The field increases with time as the charge builds up inside the Earth due to the source currents. The field reaches a steady value for times higher than the dielectric relaxation time of Earth's crust τdrt,ec=ϵ0/σec=8.9 × 10−7 s. The field reaches another steady value at the later time for times higher than the dielectric relaxation time of the atmosphere at the ground τdrt,atm=ϵ0/σ0=177 s. For the intermediate timescales 10−6−101 s the field in the highly conductive Earth is relaxed, while the weakly conducting atmosphere is not yet influencing the field. We refer to this state as an intermediate steady state due to the Earth's crust relaxation. The final state for tτdrt,atm is called the steady state due to the completion of the atmospheric relaxation. The field E0 for the case with ϵr=60, representing a higher limit of the relative permittivity of the rocks inside the Earth shows the same values as the case with ϵr=1 once both cases reach the intermediate state due to the earth's crust relaxation. This result supports an assumption about independence of the results relevant to EQLs on the relative permittivity of the earth as typical times of EQLs are longer than 10−4 s.

Figure 2a shows magnitude of the electric field for the intermediate steady state due to the Earth's crust relaxation, which corresponds to the majority of EQLs' timescales. The two-dimensional view of the charge density for this state is shown in Figure 2c. The locations of the source currents correspond to regions where the charge density is ±10−7C/m3. This value agrees with an estimate of the charge density from the source current density multiplied by the dielectric relaxation time: urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0012. Other than the charge at the sources the charge resides on the Earth's surface as is expected for a conductor. The positive surface charge ρsu is located above the dipole and forms a circle of radius of ∼10 km. The negative surface charge is located everywhere around. The weak negative charge in the atmosphere is observed above the positive surface charge, and the weak positive charge is observed above the negative surface charge. The atmospheric charges will reach steady state value on timescale tτdrt,atm and are a reason for the increase of the peak electric field E0 in Figure 2b reaching the final steady state on this timescale.

The comparison of the external energy input with the total seismic wave energy in the earthquake is shown in Figure 2d. It is seen that the energy input is increasing with time. For the intermediate steady state established due to the Earth's crust relaxation tτdrt,ec, the earthquake power is constant, and therefore the energy increases linearly with time. It is also seen that 1% of the seismic wave energy in the earthquake is used in 10−5 s. Therefore, this scenario appears to be unrealistic to produce EQLs from both the energy and the telluric currents point of view.

3.2 Case With h=−15 km and h+=0 km

In the model of Theriault et al. (2014) the charges flow from the stressed volume in the Earth's crust and they can reach the surface. In this section the upper source current location is shifted to the Earth's surface h+=0 km to examine if this scenario is feasible. The source current in the upper pole is homogeneously distributed as a surface current over a circle with of radius 1 km.

To obtain the field above the ground E0 higher than the EQLs reference threshold 105 V/m, the source current has to be set to Is=107A, 2 orders of magnitude smaller than in the previous case. The two-dimensional view of the magnitude of the electric field is shown in Figure 3a. The view corresponds to a time instant t = 1 s representing a typical timescale of EQLs. The volume where the electric field is higher than 105 V/m is enclosed by the black line, and it has a radius of 1 km and a height of few hundred meters above the Earth's surface. It encloses the upper pole of the source current. It is important to note that in Figure 2a the most of the source current can flow from one pole to the other without leaving the Earth toward the atmosphere. However, all the current flowing upward from upper positive pole leaves the Earth. Therefore, a higher fraction of source current reaches atmosphere and the lower source current Is is necessary to produce the electric field above the ground of 105 V/m than for the previous scenario. The field inside the Earth just beneath the surface at 100 km distance from the dipole is 0.6 V/m. This field is higher by several orders of magnitude than the measured values 10−6–10−4 V/m (Orihara et al., 2012; Varotsos et al., 2007).

Details are in the caption following the image
Two-dimensional view of (a) magnitude of the electric field and (c) charge density around the source currents. The volume above the Earth's surface with the electric field higher than 105 V/m is considered as the place where earthquake lights can appear. The source currents are placed at 0 and 15 km beneath the Earth's surface. The time evolution of the peak electric field E0 above the Earth's surface is shown in panel (b). The time evolution of the external work of the earthquake is shown in panel (d) and compared with the total seismic wave energy of earthquake (M = 7.5) 1016 J.

The time evolution of the peak electric field in the atmosphere E0 is shown in Figure 3b and is similar to the one observed in Figure 2b. The initial conditions correspond to a zero charge in the domain and therefore zero field. The field increases with time as charge builds up inside the Earth due to the source currents. The field reaches a steady value for times longer than the dielectric relaxation time of Earth's crust τdrt,ec=8.9 × 10−7 s. The field reaches another steady value for times longer than the dielectric relaxation time of the atmosphere near the ground τdrt,atm=177 s.

The two-dimensional view of charge density at t = 1 s is shown in Figure 3c. The positive surface charge coincides with the location of the positive source current pole. Therefore, the surface charge density can be estimated from source current surface density multiplied by the dielectric relaxation time: urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0013. The electric field above the surface is then evaluated from the formula for the field above the conductor with the surface charge density ρsu:
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0014(5)
which is in a very good agreement with the value 3.1 × 105 V/m shown in Figure 3b for the intermediate steady state.

The comparison of the external energy input with the seismic wave energy in the earthquake is shown in Figure 3d. It is observed that the electrical energy is comparable with the seismic wave energy on timescales 1–10 s. Therefore, if the source current pulses will be shorter than 10−2 s this scenario may be realistic. The decrease by 4 orders of magnitude in the energy relative to results reported in the previous section is obtained due to a decrease of the source current by 2 orders of magnitude (represented by urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0015 term in equation 8, as will be further discussed below).

3.3 Case With h=−150 m and h+=−50 m

Some EQLs observations are suggesting a smaller spatial extent of this phenomenon when compared to results in the previous sections (e.g., Derr, 1973; Heraud & Lira, 2011). Therefore, we now examine also current sources of smaller size and closer to the Earth's surface. This setup can be linked to a seismic wave propagation instead of assuming that EQLs are produced at the seismic origin. In this section the dipole's altitudes are reduced by a factor of 100 with respect to section 3. The radius of cylindrical poles in this case is 25 m, and the vertical extent of each source is 50 m.

The two-dimensional view of the magnitude of the electric field is shown in Figure 4a. The source current is Is=105A, 2 orders of magnitude smaller than in the previous case and 4 orders of magnitude smaller than in section 3. The results are shown at time instant t = 1 s representing a typical timescale of EQLs. The volume where the electric field is higher than 105 V/m is enclosed by a black line, and it has a radius of 100 m and a height of few tens of meters above the Earth's surface. The field inside the Earth just beneath the surface at a 100 km distance from the dipole is 3 × 10−6 V/m. This field is comparable with the measured values 10−6–10−4 V/m (Orihara et al., 2012; Varotsos et al., 2007).

Details are in the caption following the image
Two-dimensional view of (a) magnitude of the electric field and (c) charge density around the source currents. The volume above the Earth's surface with the electric field higher than 105 V/m is considered as the place where earthquake lights can appear. The source currents are placed at 50 and 150 m beneath the Earth's surface. The time evolution of the peak electric field E0 above the Earth's surface is shown in panel (b). The time evolution of the external work of earthquake is shown in panel (d) and compared with the total seismic wave energy of earthquake (M = 7.5) 1016 J.

The time evolution of the peak electric field in the atmosphere E0 is shown in Figure 4b and is similar to the one observed in Figure 2b. The small current in this case negligibly influences the GEC, and therefore the steady state due to atmospheric relaxation is almost equivalent to the steady state due to the Earth's crust relaxation. It also means that majority of the current in this case is directed from one pole to another pole without reaching the ionosphere and contributing to the GEC.

The two-dimensional view of the charge density for this state is shown in Figure 3c. The locations of charges inside the Earth coincide with the locations of the dipole. The positive surface charge is observed directly above the dipole and forms a circle with ∼150 m radius. The negative surface charge is located everywhere else. The positive charge at the location of the upper pole can be estimated from the source current multiplied by the dielectric relaxation time: Q+=Isτdrt,ec≃9 × 10−2C. The peak electric field E0 is then approximated from the Coulomb's law for the top positive pole:
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0016(6)
which is in very good agreement with the value 3.6 × 105 V/m shown in Figure 4b for the intermediate steady state. The Coulomb's law provides a very good estimate since the source current is so close to the surface.

The external energy input is compared with the seismic wave energy in the earthquake in Figure 4d. It is observed that the electrical energy is comparable with the total seismic wave energy on the timescale 103 s. Therefore, for the source current pulses even of duration 1 s this scenario may be realistic. The difference of 6–7 orders of magnitude in comparison with the results documented in section 3 is mostly explained by the reduction in the required source current. Constraining the charges to a sphere with radius 25 m requires a higher energy than for a sphere of radius 1.5 km leading to a 2 orders of magnitude increase. However, this increase of 2 orders is offset with decrease by 8 orders of magnitude from the reduction of the source current by 4 orders in magnitude (see further discussion related to equation 8).

This scenario may be feasible even for the higher conductivity values than σec=10−5 S/m. As will be shown later using equation 11 relating the required earthquake electrical energy with the conductivity for given required electric field, the conductivity and the earthquake electrical energy are proportional. Then, for instance, for increased conductivity σec=10−3 S/m, the required electrical energy is increased by 2 orders of magnitude, and the total seismic wave earthquake energy is reached on a timescale of ∼10 s.

4 Approximate Analytical Formulae

4.1 Formulae for Energy

In this section we develop approximate analytical formulae for the electrical energy of the source as a function of the source current and the source geometry. Then, the results are validated with the results from numerical simulations.

To estimate energy, the approximate evaluation of equation 4 is performed:
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0017(7)
where the summation is performed over both poles. The potential ϕ at the position of poles is approximated from the equation ϕ = Q/C, using the steady state charge at the dipole location Q = Isτdrt,ec=Isϵ0/σec, and capacitance C corresponds to the capacitance of the source charge region.
For a large-scale dipole from section 3 both source poles can be assumed as spheres with radius of a = 1 km. Then the capacitance of each pole is C = 4πϵ0a (Zahn, 1979, p. 178, equation (11)) leading to a formula and a value at 1 s:
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0018(8)
which agrees very well with the result 1.1 × 1019 J from the simulation shown in Figure 2d.
For a large-scale dipole from section 4 the upper pole on the Earth's surface can be considered as a disc with radius a = 1 km. Then the capacitance of the upper pole is equal to C = 8ϵ0a (Jackson, 1999, p. 136, Problem 3.3). The bottom pole is again considered as a sphere with the radius 1 km leading to an energy estimate at t = 1 s:
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0019(9)

This value is in a very good agreement with 3 × 1015 J obtained from the simulation shown in Figure 3d.

For a small-scale dipole from section 5 both source poles can be approximated as spheres with radius of a = 25 m. Then the capacitance of each pole is C = 4πϵ0a leading to a formula and a value at 1 s:
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0020(10)

This value agrees very well with 4.5 × 1012 J obtained from the simulation shown in Figure 4d.

4.2 Energy Dependence on Conductivity

The previous section considered the source current as an independent variable for evaluating the electric energy. However, in our simulation the source current is chosen to obtain the electric field above the ground to be over 105 V/m. The relation between the energy input of the earthquake and the desired electric field above the surface then provides an additional insight. Having combined equation 6 for the electric field and equation 10 for the electric energy, we obtain dependence of the electric energy as a function of desired electric field above the ground for a small-scale current dipole:
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0021(11)

Equation 11 demonstrates that for a given electric field E0 corresponding to EQLs the required earthquake energy Ws is proportional to the conductivity of the Earth σec. We note that even for other current source setups the proportionality between the required earthquake electrical energy and the conductivity remains valid. Therefore, the lower conductivity σec is more favorable to the appearance of EQLs. This provides an additional rationale for our selection of the minimum Earth's conductivity for our simulation σec=10−5 S/m.

4.3 Comparison With ULF Magnetic Measurements

In section 3 we compared the electric field measurements and the simulation results at 100 km distance from the source. In this section we focus on the magnetic field ULF measurements and compare them with an approximate analytical formulation.

The preseismic magnetic signals in the ULF range (lower than several hertz) can be observed and are reviewed in Uyeda et al. (2009, and references therein). The ULF signals have advantages over those at higher frequencies because of their large skin depth urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0022 (Balanis, 2012, section 4.3.1), where μ0 is permeability of free space and ω is the angular frequency. For a frequency of 1Hz and the conductivity value used in present work 10−5 S/m the skin depth is δ = 160 km. For higher frequencies and conductivities the skin depth is decreasing, which means that the radiation from the current pulses with shorter duration will be absorbed at shorter distances and more difficult to observe experimentally. Fraser-Smith et al. (1990) reported anomalous enhancement of the magnetic field by ∼4 nT at site located only 7 km from the epicenter. It is important to note that the time resolution of their measurements was 30 min. We make an order of magnitude estimates of the magnetic field at a distance l = 10 km from the current dipoles presented in this work based on a solution for the near field region (l/δ ≪ 1) of an infinitesimal dipole (h+hδ) (Balanis, 1997, equation (4-20d)) adjusted for an attenuation with the skin depth δ (in the considered case l/δ ≪ 1 this effect is small):
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0023
where α is the angle between the dipole and the direction toward the observer. For Is=109A and h+h=10 km the field is ∼10−2T. For Is=107A and h+h=15 km the field is ∼10−4T. For Is=105A and h+h=100 m the field is ∼10nT. In the framework of these highly simplified estimates the obtained values are in a reasonable agreement with our conclusions based on the electric field values at 100 km distance. That means that both large-scale current sources give too high values when compared with the measured ones. Only the small-scale current source gives a value comparable with the measured value. It is important to note that the large-scale current sources can still be in agreement with measured values as shown in (Bortnik et al., 2010), but with current amplitudes several orders magnitude lower than are required for the discharge initiation above the ground.

5 Discussion

5.1 Seismic Waves

Heraud and Lira (2011) reported that the luminous events correlate well with the seismic activity that occurred 150 km from the epicenter and suggested that the passage of seismic waves producing the stress on rocks can locally generate the electrical charges responsible for the local luminescence. The authors observed that the first S wave arrived 18 s before the luminescence and since then the ground acceleration was increasing until an EQL appeared. The EQLs appearance is synchronized with the maximum of the ground acceleration (meaning the highest stress levels) as is shown in Heraud and Lira (2011, Figures 10d and 10e). It is important to note that the authors observed no EQLs associated with an arrival of P waves. The precise synchronization of the passage of S wave with an increase of the electric potential on the surface inside the mine is also reported in Takeuchi et al. (2010, Figure 2). Based on the above mentioned observations, the seismic waves seem to be responsible for the local generation of the electric field and the luminous events far from the epicenter. These observations can support the results of present work, suggesting that local small size current sources require much less energy and therefore are a more probable cause of EQLs. However, it is important to emphasize that the quantitative information on a fraction of the energy transferred by the seismic waves at a distance of 150 km is not available and therefore the direct comparison is not possible.

5.2 Study of Current Pulse

In section 3 the constant source current ±Is was used. For a representation of the time-dependent source current the following function with one sine pulse is formulated:
urn:x-wiley:jgrd:media:jgrd54779:jgrd54779-math-0024(12)
where Is=109A and τ = 1 s is chosen to represent the typical EQLs' timescale. We consider the case with h=−15 km and h+=−5 km. The time evolution of the peak electric field in the atmosphere E0 is shown in Figure 2b. It is seen that at t = 0.5 s, when Is,c=Is the peak electric field is equal for both cases of constant and variable current. Under the intermediate steady state conditions the history longer than τdrt,ec does not influence the solution. If variations are on 1-s scale, which is much longer than τdrt,ec, then the intermediate steady state is established at every moment of time. The energy input of an earthquake is shown in Figure 2d. The energy input at t = 1 s is approximately half of the energy input of the constant current case. It can be concluded that for any current pulse of duration τ in the intermediate steady state regime, the order of magnitude of the energy can be estimated from a constant current case by reading the energy at the corresponding time τ.

5.3 Horizontal Size of Current Source

In this section we compare our simulation results with the results shown in Kuo et al. (2014). This scenario is characterized by a large horizontal size of the source current. Kuo et al. (2014) assume a region of 200km×450 km, which is approximated in the present work by a circular region with the radius of rKuo=170 km. The vertical locations of source currents are in the Earth's crust but are not specified in Kuo et al. (2014). For representative calculations we keep values h=−15 km and h+=−5 km. For Scur, we use the same spatial distribution as described in Kuo et al. (2014, equation (1)) for current density on the Earth's surface Jsurf, but with a = b = rKuo. We scale Scur to obtain the electric field above the center of the dipole during the intermediate steady state 3 × 104 V/m. This field corresponds to the current density 1.5nAm−2. The total upward current at the Earth's surface generated in our simulation is 82A at t = 1 s. These values match well with the lower values of the range of values provided by Kuo et al. (2014), 1–100nAm−2 and 90–9,000 A. It is important to note that the source current dipole inside the Earth necessary to generate such fields above the Earth's surface is Is=53 × 109A. The comparison of the external energy input required for EQLs based on this scenario is shown in Figure 2d. For this case, it is observed that required energy for EQLs exceeds the available earthquake energy on timescale 10−4 s, and therefore this scenario is unlikely to produce EQLs.

Kuo et al. (2015) performed a comparison of simulations and observations for the 2011 Tohoku-Oki earthquake (M = 9.0). The authors estimated the total upward current at the Earth's surface ∼2,250 A with a maximum current density 25 nA/m2. The electrical energy required for this case increases with power of 2 of the current (see equation 8), therefore ∼600 times in comparison with the reference case presented in the paragraph above. The total seismic wave energy for M = 9.0 is 1.6 × 1018 J,∼160 times higher than the reference case. The timescale when the electrical energy becomes higher than the energy of earthquake remains 10−4 s. It is important to note that the goal of Kuo et al. (2014) is to investigate the possible effects on the ionosphere by a very large current. The above estimates indicate that it is very difficult to affect the ionospheric dynamics by the lower atmospheric current source. It is also important to note that our results on energy comparison do not exclude the possibility that these large currents can exist as the energy available in the earthquake region prior to an earthquake may exceed that of the total seismic wave energy.

5.4 Conductivity Influence and Sources Beneath Sea

Section 3 describes three different scenarios of current dipoles inside the Earth with the assumed constant conductivity 10−5 S/m to represent the minimum conductivity of the Earth's crust as it provides the best case scenario for the existence of EQLs. The approximate equation 11 provides scaling of the total electrical energy as a function of conductivity for the constant electric field on the ground. It shows that the total electrical energy increases linearly with the conductivity. It is possible that a high-conductivity layer (e.g., sea) appears above the current dipole. Our simulation results show (not included here for the sake of brevity) that adding a highly conducting layer above the current dipole in section 3 decreases significantly the electric field on the surface. The high-conductivity layer keeps more current inside the Earth, and a smaller fraction of current escapes to the atmosphere, as analyzed in section 3.

5.5 Discharge Dynamics

The EQLs under the ambient electric field above 5 × 103 V/m would appear in the form of corona around trees and sharp objects (Standler & Winn, 1979). The time dynamics of nonstationary corona around a single sharp rod is discussed in Bazelyan et al. (2008, and references therein). For higher fields and sufficiently tall grounded objects an upward leader can be initiated and its characteristics are summarized in Rakov and Uman (2003, Chapter 6). It was shown that ambient potential at the tip of a grounded object has to be able to support 400kV potential drop in the streamer corona volume to be able to generate a leader (e.g., Bazelyan et al., 2008). To generate such a potential drop in the streamer corona, the potential difference between the tip of the grounded object and the ambient potential at the same altitude has to be at least twice bigger, as half of the potential drops in the streamer zone of the leader and half in free space ahead of the streamer zone (Bazelyan & Raizer, 2000, p. 69; da Silva & Pasko, 2013). For example, in the field 105 V/m, assumed in this work as a reference field, the grounded object of height 10 m would provide potential difference of 1,000 kV leading to upward leader initiation. To analyze if such leader is viable or how far it would propagate, the model developed by Becerra (2014) can be used.

6 Conclusions

The principal results of this paper are as follows:
  1. We performed quantitative investigations of various configurations of current dipoles placed inside the Earth with a purpose of explaining EQLs using an enhancement of the electric field near the Earth's surface (reference value 105 V/m) leading to conventional lightning-like gas discharges.
  2. We have also developed and tested an approximate formulation allowing effective solution of the same problems using an analytical theory.
  3. It is shown that a large-scale current dipole located in the Earth's crust with poles at 5 and 15 km beneath the Earth's surface requires during 1 s an energy of 1.1 × 1019 J significantly exceeding the total seismic wave energy 1016 J of major earthquakes of magnitude 7.5.
  4. The most likely setup to explain EQLs is found to be when the upper pole of the source current dipole is shifted close to the Earth's surface. In particular, if locations of the current dipoles are tens of meters beneath the Earth's surface, the total electrical energy required is only ∼1012 J.

Acknowledgments

This research was supported by the National Science Foundation under grants AGS-1135446 and AGS-1623780 to Penn State University. All the data used are archived in the repository doi:10.18113/D35D1G or listed in the references.