Calculation of Radon (Fluid) Precursors of Tectonic Earthquake

Igor Petrovich Dobrovolsky

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

DOI: 10.4236/oalib.1110516                

Abstract

On the Earth’s surface the continuous flow of fluids (different gases) is observed. During the preparation of an earthquake, there is the deformation of the environment which changes mass transfer conditions, and changes of this flow arise. These changes are perceived as earthquake precursor. It is the simplest to watch such variations on the flow of the radioactive gas of radon. In work calculation of this precursor in the two-phase environment is made.

Keywords

Two-Phase Medium, Gas Breath of Earth

1. Introduction

In 1912 V. I. Vernadsky at the meeting of the Russian Imperial Academy of Sciences gave the report: “About studying of gas breath of Earth”. He suggested that from the Earth interior there has to be the flow of different gases and called this phenomenon “Gas breath of Earth”. Afterwards, this assumption was confirmed in field observations. Later change of the flow of radon by preparation of tectonic earthquakes was revealed (for example, [1] ). This phenomenon received the name of the radon precursor of the tectonic earthquake. The mechanism of this precursor is almost obvious: by preparation of the tectonic earthquake there is the deformation of the porous and cracked environment (crust) which changes effective coefficient of the mass transfer.

Volume deformation has basic meaning and its distribution defines distribution of the fluid precursor. Short-term precursors because on big intervals of time there is the alignment of the flow are of special interest. It is also necessary to remember that time of half-decay of radon is about 4 days.

In this work calculation of the short-term fluid precursor for this mechanism will be made.

2. Main Notations

x, y, z is Cartesian coordinate system from the beginning in heterogeneity epicenter, [m];

t is time from the beginning of action of the short-term precursor, [s];

s is tangential stress at infinity (boundary condition), [Pa] = [kg/(m⋅s²)];

G is shear modulus of the Earth’s crust, [Pa];

α is relative increment in shear modulus in heterogeneity;

ε is the volume strain of the environment;

h is depth of hypocenter of the spherical heterogeneity, [m];

R is radius of the spherical heterogeneity, [m];

P(x, y, z, t) is the abundant porous pressure;

μ is coefficient of dynamical viscosity of fluid, [Pa⋅s];

k is permeability, [m²];

m and m₀ are the current and initial porosity of the medium;

V is velocity of fluid, [m/s];

U is velocity of solid phase, [m/s];

W is velocity of Darcy, [m/s];

β is isothermal compressibility of fluid, [1/Pa];

H is Heaviside function;

δ is delta-function;

Δ is Laplace operator with respect to spatial variables.

3. Equations of the Two-Phase Medium

In the linear approximation, in the case of slow processes, when inertial forces and compressibility of the solid phase can be neglected with respect to the compressibility of fluid, the system of equations of dynamics of a porous fluid has the form

$\mathrm{grad}P+\frac{\mu m_0}{k}\left(V-U\right)=0$ (3.1)

$\frac{1}{m_0}\frac{\partial m}{\partial t}+\beta \frac{\partial P}{\partial t}+\mathrm{div}V=0$ (3.2)

Equation (3.1) is, in essence, Darcy’s law, Equation (3.2) describes deformation of the fluid in the elementary volume of the two-phase environment.

After the transformations which are in detail stated in [2] the equation turns out

$\mathrm{div}\frac{k}{\mu m_0}\mathrm{grad}P-\beta \frac{\partial P}{\partial t}-q\frac{\partial \epsilon }{\partial t}=0$ (3.3)

In stationary mode from (3.3) we have

$\mathrm{div}\frac{k}{\mu m_0}\mathrm{grad}P=0$ (3.4)

At $k/\left(\mu m_0\right)=\mathrm{const}$ (3.3) and (3.4) give respectively the equations

$\frac{k}{\mu \beta m_0}\Delta P-\frac{\partial P}{\partial t}-\frac{q}{\beta }\frac{\partial \epsilon }{\partial t}=0$ (3.5)

$\Delta P=0$ (3.6)

4. Background Flow of Fluids

Let’s consider the background flow homogeneous and stationary. It is described by Equation (3.6) with boundary condition of free outpouring on the surface

$\frac{{\partial }^2{P}_0}{\partial z^2}=0,P_0\left(x,y,0\right)=0$ (4.1)

During the aseismatic period the flow is observed, and according to Darcy’s law is expressed as

$W_0=-\frac{k}{\mu }\mathrm{grad}{P}_0$ (4.2)

Therefore, the solution of the system (4.1) is

$P_0=\frac{\mu W_0}{k}z$ (4.3)

5. Deformation of the Environment by Preparation of the Tectonic Earthquake

From positions of mechanics of continuous medium process of preparation of the tectonic earthquake is considered as process of forming in the environment of heterogeneity of properties [3] . Two periods are allocated: long-term and short-term precursors. At the stage of long-term precursors heterogeneity arises and develops; in the period of short-term precursor it begins to collapse. At the beginning of action of short-term precursor for spherical heterogeneity at Poisson’s coefficient ν = 1/4 volume deformation has the appearance [3]

${\epsilon }_{\max }=\frac{2\alpha s{R}^3}{3G}xy\left(\frac{H\left(r_1-R\right)}{r_1^5}+\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)$ (5.1)

where $r_1=\sqrt{x^2+y^2+{\left(z-h\right)}^2}$ , $r_2=\sqrt{x^2+y^2+{\left(z+h\right)}^2}$ , $R={10}^{0.414M+1.304}$ is radius of spherical heterogeneity, M is magnitude of future earthquake.

Let’s consider the short-term precursor. Let’s assume that at the beginning of action of the short-term precursor the cumulative deformation is dumped to zero and further remains invariable

$\epsilon \left(t\right)=-{\epsilon }_{\max }H\left(t\right)$ (5.2)

Of course, also other law of behavior of deformation during action of the short-term precursor is possible.

6. The Flow Perturbed with the Earthquake Precursor

In Equation (3.5) we will denote

$a^2=\frac{k}{\mu \beta m_0}$ (6.1)

Derivative from (5.2)

$\frac{\partial \epsilon }{\partial t}=-{\epsilon }_{\max }\delta \left(t\right)$ (6.2)

Now the system from Equation (3.5) with boundary and initial conditions registers as

$\begin{array}{l}{a}^2\Delta P-\frac{\partial P}{\partial t}+\frac{q}{\beta }{\epsilon }_{\max }\delta \left(t\right)=0\\ P\left(x,y,0,t\right)=0,P\left(x,y,z,0\right)=P_0\end{array}$ (6.3)

The solution of system (6.3) is given by formula (A.2).

Important note. Thanks to the Heaviside function (5.1) suffers the gap, and at integration it leads to significant increase in time of calculations. However, such situation can be improved if to accept approximate representation

$H\left(r_1-R\right)\approx 1-\exp \left(-{\left(r_1/R\right)}^n\right)$ (6.4)

At n = 20 - 30 the error from such replacement becomes quite acceptable.

7. Merely Precursor Perturbation

To allocate such perturbation, we believe in (6.3) $P=P_0+M$ and we receive the system for M

$\begin{array}{l}{a}^2\Delta M-\frac{\partial M}{\partial t}+\frac{q}{\beta }{\epsilon }_{\max }\delta \left(t\right)=0\\ M\left(x,y,0,t\right)=0,M\left(x,y,z,0\right)=0\end{array}$ (7.1)

Then

$\begin{array}{c}M=\frac{2q\alpha s{R}^3}{3\beta G}{\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{V}{\iiint }D\left(x,y,z,t^{\prime },x^{\prime },y^{\prime },z^{\prime }\right){\epsilon }_{\mathrm{m}}\delta \left(t-t^{\prime }\right)\text{d}{x}^{\prime }\text{d}{y}^{\prime }\text{d}{z}^{\prime }\text{d}{t}^{\prime }}}\\ =\frac{2q\alpha s{R}^3}{3\beta G}{\displaystyle \underset{V}{\iiint }D\left(x,y,z,t,x^{\prime },y^{\prime },z^{\prime }\right){\epsilon }_{\mathrm{m}}\text{d}{x}^{\prime }\text{d}{y}^{\prime }\text{d}{z}^{\prime }}\end{array}$ (7.2)

In this solution dominant function is

$D\left(x,y,z,t,x^{\prime },y^{\prime },z^{\prime }\right)={\left(\frac{1}{2\sqrt{\pi a^{2}t}}\right)}^3\left({\text{e}}^{-\frac{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2+{\left(z-z^{\prime }\right)}^2}{4a^{2}t}}-{\text{e}}^{-\frac{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2+{\left(z+z^{\prime }\right)}^2}{4a^{2}t}}\right)$ (7.3)

and

${\epsilon }_{\mathrm{m}}=x^{\prime }{y}^{\prime }\left(\frac{H\left(r_1-R\right)}{r_1^5}+\frac{2}{r_2^5}-\frac{10h\left(z^{\prime }+h\right)}{r_2^7}+2R^2\left(\frac{7{\left(z^{\prime }+h\right)}^2}{r_2^9}-\frac{1}{r_2^7}\right)\right)$ (7.4)

where $r_1=\sqrt{{x^{\prime }}^2+{y^{\prime }}^2+{\left(z^{\prime }-h\right)}^2}$ , $r_2=\sqrt{{x^{\prime }}^2+{y^{\prime }}^2+{\left(z^{\prime }+h\right)}^2}$ .

Practical interest represents the fluid flow on z = 0 surface. This flow is defined by Darcy’s formula and has the appearance

$\begin{array}{c}W\left(x,y,t\right)=-\frac{k}{\mu }\frac{\partial M}{\partial z}=-\frac{2qk\alpha s{R}^3}{3\beta \mu G}{\displaystyle \underset{V}{\iiint }\frac{\partial D\left(x,y,0,t,x^{\prime },y^{\prime },z^{\prime }\right)}{\partial z}{\epsilon }_{\mathrm{m}}\text{d}{x}^{\prime }\text{d}{y}^{\prime }\text{d}{z}^{\prime }}\\ =-\frac{qk\alpha s{R}^3}{12\pi \sqrt{\pi }\beta \mu G}{\displaystyle \underset{V}{\iiint }\frac{z^{\prime }{\text{e}}^{-\frac{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2+{z^{\prime }}^2}{4a^{2}t}}}{{\left(a^{2}t\right)}^{5/2}}{\epsilon }_{\mathrm{m}}\text{d}{x}^{\prime }\text{d}{y}^{\prime }\text{d}{z}^{\prime }}\end{array}$ (7.5)

Statement and the solution of problem (7.1)-(7.5) can be executed in dimensionless quantities. We postulate

$\left(x,y,z\right)=\left(R\xi ,R\eta ,R\zeta \right),t=\frac{R^2}{a^2}\tau ,M=\frac{q}{\beta }m$ (7.6)

When transforming delta-function it is necessary to use the ratio

$\delta \left(bx\right)=\frac{1}{b}\delta \left(x\right)$ (7.7)

Then (7.1) takes the form

$\begin{array}{l}\Delta m-\frac{\partial m}{\partial \tau }+{\epsilon }_{\mathrm{m}}\left({\xi }^{\prime },{\eta }^{\prime },{\zeta }^{\prime }\right)\delta \left(\tau \right)=0\\ m\left(x,y,0,t\right)=0,m\left(x,y,z,0\right)=0\end{array}$ (7.8)

In the dimensionless form (7.2) is

$m=\frac{2\alpha s}{3G}{\displaystyle \underset{V}{\iiint }{D}^{\prime }\left(\xi ,\eta ,\zeta ,{\xi }^{\prime },{\eta }^{\prime },{\zeta }^{\prime },\tau \right){{\epsilon }^{\prime }}_{\mathrm{m}}\text{d}{\xi }^{\prime }\text{d}{\eta }^{\prime }\text{d}{\zeta }^{\prime }}$ (7.9)

where $D^{\prime }\left(\xi ,\eta ,\zeta ,{\xi }^{\prime },{\eta }^{\prime },{\zeta }^{\prime },\tau \right)={\left(\frac{1}{2\sqrt{\pi \tau }}\right)}^3\left({\text{e}}^{-\frac{{\left(\xi -{\xi }^{\prime }\right)}^2+{\left(\eta -{\eta }^{\prime }\right)}^2+{\left(\zeta -{\zeta }^{\prime }\right)}^2}{4\tau }}-{\text{e}}^{-\frac{{\left(\xi -{\xi }^{\prime }\right)}^2+{\left(\eta -{\eta }^{\prime }\right)}^2+{\left(\zeta +{\zeta }^{\prime }\right)}^2}{4\tau }}\right)$ ,

${{\epsilon }^{\prime }}_{\mathrm{m}}={\xi }^{\prime }{\eta }^{\prime }\left({\left[\frac{1}{{\rho }_1^5}\right]}_e+\frac{2}{{\rho }_2^5}-\frac{10\omega \left({\zeta }^{\prime }+\omega \right)}{{\rho }_2^7}+2\left(\frac{7{\left({\zeta }^{\prime }+\omega \right)}^2}{{\rho }_2^9}-\frac{1}{{\rho }_2^7}\right)\right)$ ,

${\rho }_1=\sqrt{{{\xi }^{\prime }}^2+{{\eta }^{\prime }}^2+{\left({\zeta }^{\prime }-\omega \right)}^2}$ , ${\rho }_2=\sqrt{{{\xi }^{\prime }}^2+{{\eta }^{\prime }}^2+{\left({\zeta }^{\prime }+\omega \right)}^2}$ , $h=R\omega$ .

Based on formula (7.5), the dimensionless speed of Darcy on the surface has the appearance

$w=-{\displaystyle \underset{V}{\iiint }\frac{{\zeta }^{\prime }{\text{e}}^{-\frac{{\left(\xi -{\xi }^{\prime }\right)}^2+{\left(\eta -{\eta }^{\prime }\right)}^2+{{\zeta }^{\prime }}^2}{4\tau }}}{{\tau }^{5/2}}{{\epsilon }^{\prime }}_{\mathrm{m}}\text{d}{\xi }^{\prime }\text{d}{\eta }^{\prime }\text{d}{\zeta }^{\prime }}$ (7.10)

Here

$W=\frac{qk\alpha s}{12\pi \sqrt{\pi }\beta \mu GR}w$ (7.11)

In Figures 1-3 lines of equal level are given in the first quadrant. On the neighboring quadrants continuation is made is antisymmetric, and, therefore, areas of decrease and increase of effect alternate. Such situation will be coordinated with empirical data.

Figure 1. The plot of function $w\left(\xi ,\eta ,\tau \right)$ при $\tau =0.01$ , $\omega =1$ in the first quadrant. Axes of coordinates are zero lines.

Figure 2. The plot of function $w\left(\xi ,\eta ,\tau \right)$ при $\tau =0.1$ , $\omega =1$ in the first quadrant. Axes of coordinates are zero lines.

Figure 3. The plot of function $w\left(\xi ,\eta ,\tau \right)$ при $\tau =1$ , $\omega =1$ in the first quadrant. Axes of coordinates are zero lines.

8. Conclusion

In formula (7.11) the radius R is in the denominator. Therefore, at increase in magnitude of future earthquake the precursor effect with other equal conditions will be less. Empirically it is almost impossible to check this conclusion as it is difficult to meet “other equal conditions”. On the other hand, this conclusion is not unexpected. For the earthquake of bigger magnitude, the same deformation extends to the big area and becomes more flat.

The main conclusion can be formulated so: the chosen technique can be applied to calculations of fluid precursors of the tectonic earthquake.

Appendix: Parabolic Equation

The first boundary value problem for the parabolic equation with constant coefficients and its solution for half-spaces are defined by the following formulas

$\begin{array}{l}{a}^2\Delta u-\frac{\partial u}{\partial t}+f=0\\ u\left(x,y,0,t\right)=\Phi \left(x,y,t\right),u\left(x,y,z,0\right)=\psi \left(x,y,z\right)\end{array}$ (A.1)

and the solution is expressed by means of dominant function

$u\left(x,y,z,t\right)=a^2{\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{S}{\iint }\Phi \frac{\partial D}{\partial z^{\prime }}\text{d}{s}^{\prime }\text{d}{t}^{\prime }}+{\displaystyle \underset{V}{\iiint }\psi D\text{d}{v}^{\prime }}+{\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{V}{\iiint }Df\text{d}{v}^{\prime }\text{d}{t}^{\prime }}}}$ (A.2)

where S and V are the surface and volume of the half-space respectively, integration on time is made in the form of convolution.

Dominant function for half-spaces $z\ge 0$

$D\left(x,y,z,t,x^{\prime },y^{\prime },z^{\prime }\right)={\left(\frac{1}{2\sqrt{\pi a^{2}t}}\right)}^3\left({\text{e}}^{-\frac{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2+{\left(z-z^{\prime }\right)}^2}{4a^{2}t}}-{\text{e}}^{-\frac{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2+{\left(z+z^{\prime }\right)}^2}{4a^{2}t}}\right)$ (A.3)

In the analysis of this solution, it is useful to mean ratio

$\underset{t\to 0}{\lim }\frac{1}{2\sqrt{\pi t}}{\text{e}}^{-\frac{x^2}{4t}}=\delta \left(x\right)$ (A.4)

These data can be found in any, rather full textbook on the equations of mathematical physics or in the reference book on mathematics.

Conflicts of Interest

The author declares no conflicts of interest.

References

[1]	Rikitake, T. (1976) Earthquake Prediction. Elsevier Scientific Publishing Company, Amsterdam.
[2]	Dobrovolsky, I.P. (2014) Hydrodynamical Phenomena at Preparation of Earthquake. Open Access Library Journal, 1, e896. http://dx.doi.org/10.4236/oalib.1100896
[3]	Dobrovolsky, I.P. (2022) Heterogeneity in the Elastic Half-Space (Deformations at Preparation of the Tectonic Earthquake). Open Access Library Journal, 9, e8714. https://doi.org/10.4236/oalib.1108714