The Precursor of Tectonic Earthquake in Telluric Current

Igor Petrovich Dobrovolsky

Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russia.

DOI: 10.4236/oalib.1103488                

Abstract

The calculation of the precursor of the tectonic earthquake in the telluric current is carried out on the basis of the main equations of the telluric current and the theory of preparation of the tectonic earthquake. Results do not contradict data of field observations. Methodical significance of work is especially important.

Keywords

The Solution of Equation of Telluric Current, The Deformation at Earthquake’s Preparation

1. Introduction

Field observations of electric precursors of a tectonic earthquake are part of the general observations over telluric currents and are carried out long ago. However, there are no well-founded theoretical works on this subject still. To calculate such precursor two components are necessary: The equations of telluric current and effective theory of a tectonic earthquake’s preparation. The equations of telluric current are consequence of Maxwell’s equations and the corresponding material equations. These data are known and entered tens of textbooks and reference books. The theory of preparation was constructed later [1] [2] and is a little used in the analysis of processes of earthquake preparation.

The submitted article liquidates this gap. Article purpose is to receive real estimates of size of a precursor in telluric current. Calculations are made in SI.

2. Main Equation of Telluric Current

Main equation of telluric current is based on the continuity equation which is a consequence of the equations of Maxwell, and has form

$\text{div}j+\frac{\partial \rho }{\partial t}=0$ (2.1)

where j is current density, A/m²; $\rho$ is volume density of electric charges, C/m³; t is time.

The Ohm’s law in which electric field strength is expressed through potential is used

$j=\sigma E\text{or}j=-\sigma \text{grad}V$ (2.2)

where E―electric field strength, V/m; $\sigma$ ―specific conductivity, S/m; V―potential, V.

Bringing the second Equations ((2.2) in (2.1)), we receive the known equation for time-constant telluric current

$\sigma {\nabla }^{2}V+\mathrm{grad}\sigma \cdot \mathrm{grad}V=0$ (2.3)

where ${\nabla }^2$ is Laplace’s operator.

3. Change of Conductivity at Earthquake Preparation

Deformation of the medium causes conductivity change. We will assume that change of conductivity is proportional to volume deformation

$\Delta \sigma =-\bar{\sigma }\gamma$ (3.1)

where $\gamma$ is volume deformation, $\bar{\sigma }$ is proportionality coefficient, S/m.

The sign “minus” in (3.1) reflects the following assumption: the volume of the element of the medium increases at positive deformation and conductivity decreases. In this case $\bar{\sigma }$ is positive.

Process of preparation of a tectonic earthquake consists in formation of heterogeneity in a bowels of the Earth [2] . It is always followed by deformation of the medium. Let’s consider the problem in the Cartesian axials (x, y, z) in half-space z ≥ 0. Let’s accept heterogeneity in the form of the sphere. In these conditions volume deformation of the medium is expressed by the formula [2]

$\gamma =\frac{\alpha \tau R^3\left(1-2\nu \right)}{\mu \left(1-\nu \right)}xy\left[\left(\frac{1}{r_1^5}\right){\delta }_e+\frac{3-4\nu }{r_2^5}-\frac{10H\left(z+H\right)}{r_2^7}+2R^2\left(\frac{7{\left(z+H\right)}^2}{r_2^9}-\frac{1}{r_2^7}\right)\right]$ (3.2)

where R is radius of spherical heterogeneity, μ is shear modulus out of heterogeneity, $\alpha$ is relative change of the shear modulus in heterogeneity, τ are shearing stresses on infinity, H is depth of the center of heterogeneity, ν is Poisson’s coefficient; δ_e is characteristic function of area out of heterogeneity (δ_e = 1 out of heterogeneity and δ_e = 0 in heterogeneity), $r_{1,2}=\sqrt{x^2+y^2+{\left(z\mp H\right)}^2}$ ; shearing stresses τ are directed on axes (x, y).

Accepting $\nu =1/4$ , we receive from (3.1) and (3.2)

$\begin{array}{l}\Delta \sigma =-\frac{2\alpha \tau R^3\bar{\sigma }}{3\mu }\epsilon \\ =-\frac{2\alpha \tau R^3\bar{\sigma }}{3\mu }xy\left[\left(\frac{1}{r_1^5}\right){\delta }_e+\frac{2}{r_2^5}-\frac{10H\left(z+H\right)}{r_2^7}+2R^2\left(\frac{7{\left(z+H\right)}^2}{r_2^9}-\frac{1}{r_2^7}\right)\right]\end{array}$ (3.3)

where designation $\epsilon$ is obvious.

Radius of the sphere accepts the maximum value before the earthquake

$R_{\max }={10}^{0.414M+1.304}m$ (3.4)

where M is the magnitude of future earthquake. Obviously $H\ge R_{\max }$ .

4. Statement and Solution of the Problem

We consider in half-space z ≥ 0 two problems. In the first problem the telluric current is described prior to earthquake preparation. In the second problem initial conductivity is broken by process of preparation of the earthquake and, therefore, there is the redistribution of the telluric current. The difference of solutions of these problems is the earthquake harbinger in the telluric current. Boundary conditions in both problems are identical.

We assume that in the first problem conductivity is constant $\sigma \left(x,y,z\right)={\sigma }_0=\mathrm{const}$ and in the medium homogeneous current flows $j_0\left(x,y,z\right)=\left\{j_0,0,0\right\}=\mathrm{const}$ . On the half-space surface current does not proceed in the atmosphere. Proceeding from (2.3), such problem is described by system of equations for potential $\phi \left(x,y,z\right)$

$\begin{array}{l}{\nabla }^2\phi =0\\ j\left(\pm \infty ,\pm \infty ,z\right)=j_0=\left\{j_0,0,0\right\}=\mathrm{const}\\ \frac{\partial \phi \left(x,y,0\right)}{\partial z}=0\end{array}$ (4.1)

The system (4.1) has the solution

$\phi =-\frac{j_0}{{\sigma }_0}x$ (4.2)

At earthquake preparation conductivity changes at the value determined by the formula (3.3)

${\sigma }_e={\sigma }_0+\Delta \sigma$ (4.3)

Thus, potential $\psi$ at preparation of the earthquake is defined by system

$\begin{array}{l}{\sigma }_e{\nabla }^2\psi +\mathrm{grad}\Delta \sigma \cdot \mathrm{grad}\psi =0\\ j\left(\pm \infty ,\pm \infty ,z\right)=j_0=\left\{j_0,0,0\right\}=\mathrm{const}\\ \frac{\partial \psi \left(x,y,0\right)}{\partial z}=0\end{array}$ (4.4)

If to designate $q=\psi -\phi$ , then for q we have system

$\begin{array}{l}{\sigma }_e{\nabla }^{2}q+\mathrm{grad}\Delta \sigma \cdot \mathrm{grad}\phi +\mathrm{grad}\Delta \sigma \cdot \mathrm{grad}q=0\\ j\left(\pm \infty ,\pm \infty ,z\right)=0\\ \frac{\partial q\left(x,y,0\right)}{\partial z}=0\end{array}$ (4.5)

Let’s consider also the task for half-space z ≥ 0

$\begin{array}{l}{\nabla }^{2}G=\delta \left(x-\xi ,y-\eta ,z-\zeta \right)\\ \frac{\partial G\left(x,y,0,\xi ,\eta ,\zeta \right)}{\partial z}=0\end{array}$ (4.6)

where δ is Dirac’s function.

The solution of such task is source function

$\begin{array}{l}G\left(x,y,z,\xi ,\eta ,\zeta \right)\\ =\frac{F\left(x,y,z,\xi ,\eta ,\zeta \right)}{4\text{π}}\\ \text{=}\frac{1}{4\text{π}}\left(\frac{1}{\sqrt{{\left(x-\xi \right)}^2+{\left(y-\eta \right)}^2+{\left(z-\zeta \right)}^2}}+\frac{1}{\sqrt{{\left(x-\xi \right)}^2+{\left(y-\eta \right)}^2+{\left(z+\zeta \right)}^2}}\right)\end{array}$ (4.7)

Let’s bring (4.2) in the first Equation (4.5)

${\nabla }^{2}q+\frac{1}{{\sigma }_e}\mathrm{grad}\Delta \sigma \cdot \mathrm{grad}q=\frac{j_0}{{\sigma }_0{\sigma }_e}\frac{\partial \left(\Delta \sigma \right)}{\partial x}$ (4.8)

Thus, the system (4.5), (4.8) in view of (4.6), (4.7) is equivalent to the integro-differential equation of the II kind

$\begin{array}{l}q\left(x,y,z\right)+{\displaystyle \underset{z\ge 0}{\iiint }\frac{G\left(x,y,z,\xi ,\eta ,\zeta \right)}{{\sigma }_0+\Delta \sigma }\mathrm{grad}\Delta \sigma \cdot \mathrm{grad}qdv\left(\xi ,\eta ,\zeta \right)}\\ =\frac{j_0}{{\sigma }_0}{\displaystyle \underset{z\ge 0}{\iiint }\frac{G\left(x,y,z,\xi ,\eta ,\zeta \right)}{{\sigma }_0+\Delta \sigma }\frac{\partial \left(\Delta \sigma \left(\xi ,\eta ,\zeta \right)\right)}{\partial \xi }\text{d}v\left(\xi ,\eta ,\zeta \right)}\end{array}$ (4.9)

Deformations by preparation of the earthquake have small value. According to field observations [3] the maximum value of deformation has 10⁻⁴ in future focal area and decreases with distance from epicenter. Therefore $\left|\Delta \sigma \right|\ll {\sigma }_0$ . Neglecting $\Delta \sigma$ in comparison with ${\sigma }_0$ , from (4.9) we will receive simpler equation

$\begin{array}{l}q\left(x,y,z\right)+\frac{1}{{\sigma }_0}{\displaystyle \underset{z\ge 0}{\iiint }G\left(x,y,z,\xi ,\eta ,\zeta \right)\mathrm{grad}\Delta \sigma \cdot \mathrm{grad}q\text{d}v\left(\xi ,\eta ,\zeta \right)}\\ =\frac{j_0}{{\sigma }_0^2}{\displaystyle \underset{z\ge 0}{\iiint }G\left(x,y,z,\xi ,\eta ,\zeta \right)\frac{\partial \left(\Delta \sigma \left(\xi ,\eta ,\zeta \right)\right)}{\partial \xi }\text{d}v\left(\xi ,\eta ,\zeta \right)}\end{array}$ (4.10)

The Equations ((4.9) and (4.10)) are difficult and it is impossible to receive their exact solution. As we are interested in estimate of function q and there is the small parameter in (4.9) and (4.10) $\left|\Delta \sigma /{\sigma }_0\right|\ll 1$ , we will use the null approximation in the solution of the Equation (4.10)

$q\left(x,y,z\right)=\frac{j_0}{{\sigma }_0^2}{\displaystyle \underset{z\ge 0}{\iiint }G\left(x,y,z,\xi ,\eta ,\zeta \right)\frac{\partial \left(\Delta \sigma \left(\xi ,\eta ,\zeta \right)\right)}{\partial \xi }\text{d}v\left(\xi ,\eta ,\zeta \right)}$ (4.11)

Now the task consists in calculation of integral (4.11). If in (4.11) to make the integration by parts on $\xi$ , then we will receive

$q\left(x,y,z\right)=-\frac{j_0}{4\text{π}{\sigma }_0^2}{\displaystyle \underset{z\ge 0}{\iiint }\frac{\partial F\left(x,y,z,\xi ,\eta ,\zeta \right)}{\partial \xi }\Delta \sigma \left(\xi ,\eta ,\zeta \right)\text{d}v\left(\xi ,\eta ,\zeta \right)}$ (4.12)

We bring in (4.12) value Δσ from (3.3)

$q\left(x,y,z\right)=\frac{\alpha \tau \bar{\sigma }{j}_0}{6\text{π}\mu {\sigma }_0^2}{\displaystyle \underset{z\ge 0}{\iiint }\frac{\partial F\left(x,y,z,\xi ,\eta ,\zeta \right)}{\partial \xi }{R}^3\epsilon \left(\xi ,\eta ,\zeta \right)\text{d}v\left(\xi ,\eta ,\zeta \right)}$ (4.13)

We enter into (4.13) dimensionless coordinates $\left(u,v,w\right)=\left(x/R,y/R,z/R\right)$ and we make the change of integration variables $\left(\xi ,\eta ,\zeta \right)=\left(\alpha ,\beta ,\chi \right)=\left(\xi /R,\eta /R,\zeta /R\right)$ . It means: we enter new integration variables and designate them old letters. As data on the Earth’s surface have practical value, we suppose w = 0 (z = 0). As a result of (4.13) receives the form

$q\left(u,v,0\right)=\frac{\alpha \tau \bar{\sigma }{j}_{0}R}{3\text{π}\mu {\sigma }_0^2}{\displaystyle \underset{z\ge 0}{\iiint }f\left(u,v,\xi ,\eta ,\zeta \right)\epsilon \left(\xi ,\eta ,\zeta \right)\text{d}v\left(\xi ,\eta ,\zeta \right)}$ (4.14)

where

$\begin{array}{l}f\left(u,v,\xi ,\eta ,\zeta \right)=\frac{u-\xi }{{\left({\left(\xi -u\right)}^2+{\left(\eta -v\right)}^2+{\zeta }^2\right)}^{3/2}}\\ \epsilon =\xi \eta \left(\frac{{\delta }_e}{{\rho }_1^5}+\frac{2}{{\rho }_2^5}-\frac{10h\left(\zeta -h\right)+2}{{\rho }_2^7}+\frac{14{\left(\zeta +h\right)}^2}{{\rho }_2^9}\right)\\ {\rho }_{1,2}=\sqrt{{\xi }^2+{\eta }^2+{\left(\zeta \mp h\right)}^2},h=H/R\ge 1\end{array}$

The integral (4.14) has singularity in the point $\left(\xi =u,\eta =v,\zeta =0\right)$ . Therefore, we do change of the integration variables $\left(\alpha ,\beta ,\zeta \right)=\left(\xi -u,\eta -v,\zeta \right)$ . In new variables we pass to spherical system $\alpha =\rho \sin \theta \cos \phi$ , $\beta =\rho \sin \theta \sin \phi$ , $\zeta =\rho \cos \theta$ , in which the volume element is $\text{d}v={\rho }^2\sin \theta \text{d}\phi \text{d}\theta \text{d}\rho$ . Also we will consider surface earthquakes at h = 1 when the effect is maximum. For fast approximate calculation of integral in (4.14) there is one more obstacle: characteristic function ${\delta }_e$ does the integrand discontinuous. It is possible to bypass this trouble in two ways. First, it is possible to make calculation in two areas, but it increases bulkiness of calculations, and secondly, to approximate function ${\delta }_e$ by smooth function. Obviosly

$1-\exp \left(-{\rho }_1^n\right)\to {\delta }_eприn\to \infty$ (4.15)

We will use expression (4.15) at n = 30.

As the result the integral from (4.14) receives the form

$\omega ={\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \underset{0}{\overset{\text{π}/2}{\int }}{\displaystyle \underset{-\text{π}}{\overset{\text{π}}{\int }}\Phi \cdot \chi \cdot {\rm E}\text{d}\phi \text{d}\theta \text{d}\rho }}}$ (4.16)

where

$\begin{array}{l}\Phi =-\cos \phi {\left(\sin \theta \right)}^2\\ \chi =\xi \eta =0.5{\rho }^2{\left(\sin \theta \right)}^2\sin 2\phi +\rho \sin \theta \left(v\cos \phi +u\sin \phi \right)+uv\\ {\rm E}=\frac{1-\exp \left(-{\rho }_1^{30}\right)}{{\rho }_1^5}+\frac{2}{{\rho }_1^5}-\frac{10\rho \cos \theta +8}{{\rho }_1^7}+\frac{14{\left(\rho \cos \theta +1\right)}^2}{{\rho }_1^9}\\ {\rho }_{1,2}=\sqrt{u^2+v^2+1+{\rho }^2+2\rho \left(u\cos \phi \sin \theta +v\sin \theta \sin \phi \mp \cos \theta \right)}\end{array}$

Expression (4.16) enters the new dimensionless function and we have

$q\left(u,v,0\right)=\frac{\alpha \tau \bar{\sigma }{j}_{0}R}{3\text{π}\mu {\sigma }_0^2}\omega \left(u,v\right)$ (4.17)

Isolines of function ω are given in Figure 1 and they give the presentation of changes of potential in connection with earthquake preparation.

Figure 1. Isolines of dimensionless potential ω in the first quatrant.

5. Results and Discussions

By means of formulas (4.16), (4.17) and drawing of isolines it is possible to carry out numerical estimates of effect. For this purpose at first it is necessary to transform the formula (4.17) so that to make comparison in comparable units for different earthquakes. Let we measured $\Delta \omega$ on dimensionless base $\Delta \rho$ . If to separate $\Delta q$ on $R\Delta \rho$ and to increase on 1000, then we will receive $\Delta q$ in V/km

$\Delta q=\frac{1000\alpha \tau \bar{\sigma }{j}_0}{3\text{π}\mu {\sigma }_0^2\Delta \rho }\Delta \omega \text{V}/\text{km}$ (5.1)

Conclusion follows from the formula (5.1): the value of effect does not depend on the magnitude of future earthquake with other things being equal

Now it is necessary to define coefficient $\bar{\sigma }$ . In [4] it is suggested that change of conductivity at preparation of the earthquake does not exceed 20%. In [3] it is told that on empirical data the maximum deformation at preparation of the earthquake does not exceed ${10}^{-4}$ and it matches calculations for the formula (3.3). Considering these data we suppose

$\bar{\sigma }=1000k{\sigma }_0$ (5.2)

On this formula at the maximum deformation and k = 1 we receive change of conductivity for 10%. The value of k should not exceed 3. As result from (5.1) and (5.2) we have

$\Delta q=\frac{{10}^6\alpha \tau k{j}_0}{3\text{π}\mu {\sigma }_0\Delta \rho }\Delta \omega \text{V}/\text{km}$ (5.3)

Follows from Figure 1 that the maximum value $\Delta \omega =0.14$ is reached at $\Delta \rho =0.1$ . We accept also: $\alpha =0.1$ , $\tau ={10}^8\text{Pa}$ , $\mu =2\times {10}^{10}\text{Pa}$ , k = 1, $j_0=2\times {10}^{-6}\text{A}/{\text{m}}^{\text{2}}$ , ${\sigma }_0=0.01\text{S}/\text{m}$ . As result we receive Δq = 15 mV/km. Let’s remind that conductivity of crust changes over a wide range: 10³ - 10^-7 S/m. The reader can carry out calculation with other set of parameters. We consider that detection of precursors of the tectonic earthquake in the telluric current is possible.

Acknowledgements

The author expresses sincere gratitude of E. Yu. Sokolova who reported to the author data on researches in the field of telluric currents.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Dobrovolsky, I.P., Zubkov, S.I. and Miachkin, V.I. (1979) Estimation of the Sizes of Earthquake Preparation Zones. Pageoph, 117, 1025-1044. https://doi.org/10.1007/BF00876083
[2]	Dobrovolsky, I.P. (2009) The Mathematical Theory of Preparation and Prediction of Tectonic Earthquake. Moscow. FIZMATLIT, 240 p. (in Russian).
[3]	Kasahara, K. (1981) Earthquake Mechanics. Cambrige Univ. Press, 253 p.
[4]	Rikitake, T. (1976) Earthquake Prediction. Elsevier Scientific Publishing Company. Amsterdam, 381 p.