*Editorial* **Special Issue on Comprehensive Research in Earthquake Forecasting and Seismic Hazard Assessment**

**Alexey Zavyalov 1,\* and Eleftheria Papadimitriou 2,\***


Dear Colleagues,

Despite some success, the issue of earthquake forecasting has yet to be resolved. There are occasional discussions within the scientific community about the principal feasibility of earthquake forecasting, particularly in the short-term aspect. However, the bulk of these discussions were set in the Resolution of the General Assembly of the International Association of Seismology and Physics of the Earth's Interior (IASPEI) in 2009 in Cape Town: "Resolution 4: Earthquake Forecasting and Predictability Studies—IASPEI RECOG-NIZING the opportunities provided by recent developments in earthquake science and technology RECOMMENDS that research on forecasting and predictability of earthquakes, and the validation and comparative testing of prediction methods be supported".

However, it is not sufficient to precisely predict a future strong earthquake. It is necessary to make a correct, scientifically based assessment of the level of seismic hazard and the intensity of seismic shocks to be expected in a particular region, city and settlement. What should the administration of a megapolis do when it receives information about the likelihood of a strong earthquake? The problems of earthquake forecasting and seismic hazard assessment are, therefore, closely related to the problems of high-quality anti-seismic constructions.

More than 13 years have passed since the adoption of the IASPEI Resolution. New earthquakes have occurred. Their study increased our knowledge regarding the physics of the seismic process, the physics of earthquake preparation processes and the search for earthquake precursors. The new data obtained became the basis for the development of new models of the behaviour of the ground under the influence of seismic waves and provided initial information for the development and parameterization of earthquake occurrence zone models and ground motion prediction equations.

More than one and a half years have passed since the announcement of the Special Issue "Comprehensive Research in Earthquake Forecasting and Seismic Hazard Assessment" in the MDPI Journal of *Applied Sciences*. We invited representatives of the seismological community to present their results on these topics, to show the current view of the state of the problem, what has been achieved in the field of earthquake forecasting and seismic hazard assessment, what needs to be done next and in which direction to move forward. We expected to discuss the results and directions of further research on the physics of the seismic process—from experiments under laboratory conditions to rock bursts in mines and earthquakes in seismically active regions at the stage of preparation for strong earthquakes.

As a result, 14 articles were published in the Special Issue, with authors representing different thematic areas and working in different institutions and organisations in Russia, Greece, Italy, Colombia, New Zealand, China, Argentina and Japan. The total number of authors was around 50. Thus, we managed to attract a sufficiently wide range of representatives of the scientific geophysical community to participate in this Special Issue. In this sense, our hopes and assumptions were fulfilled. In addition, this Special Issue is

**Citation:** Zavyalov, A.;

Papadimitriou, E. Special Issue on Comprehensive Research in Earthquake Forecasting and Seismic Hazard Assessment. *Appl. Sci.* **2023**, *13*, 11564. https://doi.org/10.3390/ app132011564

Received: 18 October 2023 Accepted: 19 October 2023 Published: 23 October 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

fully in line with Resolution 3 "Sharing Geophysical Data across Borders" adopted at the 28th IUGG General Assembly in Berlin on 18 July 2023.

All published articles can be roughly divided into three unequal groups in terms of the number of articles presented. The first group includes theoretical and methodological articles [1,2]. The second group includes articles confirming one or another model of seismicity behaviour in anticipation of a strong earthquake [3–5]. Finally, the third and most numerous group of articles consists of those analysing the results of long-term observations of the behaviour of various geophysical fields (seismic noise [6], seismicity [7–9], magnetotelluric field [10,11], deformation field [12], infrared radiation [13], vertical electric field in the atmosphere [14]) before strong earthquakes. We are confident that each of these articles will find an interested reader, and the whole collection will deserve the attention of representatives of the scientific community dealing with the problem of earthquake forecasting and the search for their precursors.

**Acknowledgments:** We are grateful to all the contributors who made this Special Issue, "Comprehensive Research in Earthquake Forecasting and Seismic Hazard Assessment", a success. We thank and congratulate all the authors for submitting their papers. Our sincere gratitude is also extended to all the reviewers for their efforts and time spent in helping the authors to improve their papers. We would like to express our gratitude to the Editorial Team of *Applied Sciences* for their effective and tireless editorial support for the success of this Special Issue.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **On the Omori Law in the Physics of Earthquakes**

**Alexey Zavyalov 1,\*, Oleg Zotov 1,2, Anatol Guglielmi <sup>1</sup> and Boris Klain <sup>2</sup>**


**Abstract:** This paper proposes phenomenological equations that describe various aspects of aftershock evolution: elementary master equation, logistic equation, stochastic equation, and nonlinear diffusion equation. The elementary master equation is a first-order differential equation with a quadratic term. It is completely equivalent to Omori's law. The equation allows us to introduce the idea of proper time of earthquake source "cooling down" after the main shock. Using the elementary master equation, one can pose and solve an inverse problem, the purpose of which is to measure the deactivation coefficient of an earthquake source. It has been found for the first time that the deactivation coefficient decreases with increasing magnitude of the main shock. The logistic equation is used to construct a phase portrait of a dynamical system simulating the evolution of aftershocks. The stochastic equation can be used to model fluctuation phenomena, and the nonlinear diffusion equation provides a framework for understanding the spatiotemporal distribution of aftershocks. Earthquake triads, which are a natural trinity of foreshocks, main shock, and aftershocks, are considered. Examples of the classical triad, the mirror triad, the symmetrical triad, as well as the Grande Terremoto Solitario, which can be considered as an anomalous symmetrical triad, are given. Prospects for further development of the phenomenology of earthquakes are outlined.

**Keywords:** earthquake; source deactivation; logistic equation; nonlinear diffusion equation; Omori epoch; round-the-world echo; mirror triad

#### **1. Introduction**

Omori Law [1] describes the evolution of the aftershocks of a strong earthquake. Established at the end of the century before last, the law is characterized by the beauty of its form, quite definite clarity, as a result of which it still attracts considerable attention from the geophysical community (e.g., see [2–4]).

Initially, Omori law, which can be called hyperbolic, was formulated as follows:

$$m(t) = k/(c+t).\tag{1}$$

Here *n* is the frequency of aftershocks [1]. Formula (1) is one-parameter since the parameter *c* is free and is completely determined by an arbitrary choice of the time origin. In contrast to the aspirations of Hirano [5] and Utsu [6], who introduced a two-parameter modification of Formula (1) into widespread use, we came to the conclusion that it is reasonable to put the differential equation of evolution into the basis of the phenomenological theory of aftershocks [7–9]. Guided by this consideration, in recent years we have accumulated considerable experience in the study of the evolution of aftershocks.

The purpose of this paper is to summarize our results. The main attention is focused on the phenomenological theory of aftershocks. The paper indicates successful examples of the use of theory in the analysis of experimental material. We also paid some attention to the presentation of our position on controversial issues of a methodological nature. For

**Citation:** Zavyalov, A.; Zotov, O.; Guglielmi, A.; Klain, B. On the Omori Law in the Physics of Earthquakes. *Appl. Sci.* **2022**, *12*, 9965. https:// doi.org/10.3390/app12199965

Academic Editor: Nicholas Vassiliou Sarlis

Received: 14 August 2022 Accepted: 30 September 2022 Published: 4 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

example, in the literature there are erroneous statements about the deep physical content of the parameter *c* in Formula (1).

Another misconception is associated with the idea that the two-parameter Hirano– Utsu formula is preferable to the one-parameter Omori formula. As an argument, it is argued that the presence of two parameters facilitates the approximation of observation data. On the contrary, we consider the presence of only one physical parameter in Formula (1) as an important sign of the fundamental nature of the Omori law.

#### **2. Elementary Master Equation**

The differential approach to modeling aftershocks opens up a wide scope for searches. From the richest set of differential equations available here, we take the simplest implementation of our idea, namely the truncated Bernoulli equation:

$$d n/dt + \sigma n^2 = 0,\tag{2}$$

Here *σ* is the so-called deactivation coefficient of the earthquake source, "cooling down" after the main shock [7–9]. Elementary master Equation (2) is useful in making the law of evolution simpler and easier to understand. It expresses the essence of the Omori hyperbolic law (1) that everyone understands. Moreover, it will serve as an initial basis for interesting generalizations (see below Sections 3–5).

Equation (2) contains only one phenomenological parameter *σ*. It is easy to make sure that both formulations of the law, (1) and (2), are completely equivalent to each other for *σ* = *const*. However, firstly, in contrast to (1), Formula (2) makes it possible to take into account the nonstationarity of the geological medium in the source, which undergoes a complex relaxation process after the discontinuity has formed during the main shock. The second advantage of Formula (2) is no less important. We can seek and find natural generalizations of the differential law of evolution of aftershocks, which opens up new, sometimes unexpected, approaches to processing and analyzing experimental data.

Concluding this section of the paper, let us show how easy it is to take into account nonstationarity when formulating the Omori law in the form of the evolution Equation (2). For this, it is sufficient to assume that the deactivation factor depends on time. Let us rewrite Omori's law in the most compact form:

$$dn/d\tau + n^2 = 0,\tag{3}$$

where *τ* = *t* <sup>0</sup> *σ*(*t* - )*dt*- . The general solution to Equation (3) is

$$n(\pi) = n\_0/(1 + n\_0 \pi),\tag{4}$$

It is seen that solution (4) retains the hyperbolic structure of the law, which was originally established thanks to Omori's discernment. The difference between (4) and (1) is only that time in the source, figuratively speaking, flows unevenly. For *σ* = *const*, (4) coincides with (1) up to notation. Thus, Equation (2) and its solution (4), in a certain sense, complete Omori's plan, as well as what Hirano and Utsu were striving for in their attempt to improve the law using a two-parameter modification of Formula (1).

#### **3. Logistic Equation**

Faraoni [10] considered the possibility of representing the Omori law, written in the form (2), as the Euler–Lagrange equation. The Lagrangian formulation of Omori law is interesting in many ways. In particular, it provides a basis for searching for possible generalizations of the law [11]. But the search for a suitable Lagrangian is not the only way to derive the evolution equation. One can proceed, for example, from a fairly general integro-differential equation:

$$\frac{dn}{dt} = \int\_0^\infty K(t - t') \cdot F[n(t')] \cdot dt'.\tag{5}$$

If we put *<sup>F</sup>*(*n*) = −*σn*2, then, when choosing the trivial kernel *<sup>K</sup>*(*<sup>t</sup>* − *<sup>t</sup>* - ) = *δ*(*t* − *t* - ) from (5) follows the Omori law in the form (2).

Derivation of (2) from (5) is useful in the sense that a natural generalization of the Omori law is suggested to us. The need for generalization is dictated by the following consideration. It follows from (2) that lim*n*(*t*) = 0 for *t* → ∞. Meanwhile, experience shows that the flow of aftershocks ends with a transition to a certain background seismicity of the source. It is desirable for us to take this circumstance into account by using minimal changes in the form of the classical Omori law. It turns out that for this it is enough to take into account the linear term in the formula *<sup>F</sup>*(*n*): *<sup>F</sup>*(*n*) = *<sup>γ</sup><sup>n</sup>* − *<sup>σ</sup>n*2. Here *<sup>γ</sup>* is the second phenomenological parameter of our theory. As a result, we get the master equation in the following form:

$$\frac{dn}{dt} = n(\gamma - \sigma n). \tag{6}$$

This is the logistic equation of Verhulst [12], well known in biology, chemistry, and sociology. It turns out to be useful in the physics of earthquakes [9,11].

We divide the family of solutions to logistic Equation (6) into two classes. The first class includes growing, and the second, falling functions of time. The separation principle is easiest to show in the phase portrait shown in Figure 1, where we have used the dimensionless quantities

$$X = \frac{n}{n\_{\text{max}}}, \ P = \frac{n\_{\infty}}{\gamma n\_{\text{max}}} \cdot \frac{dX}{dT}, \ T = \gamma t \tag{7}$$

instead of the original quantities. Here *n*<sup>∞</sup> = *γ*/*σ*.

**Figure 1.** Phase portrait on the phase plane of Equation (6). The red, green, and blue phase trajectories are plotted at *X*∞ = 0, 0.5 and 1, respectively (see text).

Faraoni [10] proposed the introduction of the phase plane of the dynamic system simulating the evolution of aftershocks according to the Omori law (2). The corresponding phase portrait is shown in Figure 1 with the red line. The point (0,0) corresponds to the equilibrium state. The representative point moves from bottom to top along the phase trajectory with deceleration. The portrait consists of one phase trajectory that starts at point (1,−1) and ends at point (0,0).

However, we are interested in the family of phase trajectories for Equation (6), constructed for different values of the parameter *X*<sup>∞</sup> = *n*∞/*nmax*. The red, green, and blue trajectories in Figure 1 are plotted with *X*∞ values of 0, 0.5, and 1, respectively. Equilibrium point (0,0) is stable at *X*<sup>∞</sup> = 0 and unstable at *X*<sup>∞</sup> > 0. Equilibrium point (*X*∞, 0) is stable at any values of parameter *X*<sup>∞</sup> > 0. It can be shown that the velocity of motion of the imaging point along the phase trajectory asymptotically tends to zero with approaching (*X*∞, 0). The segment of the trajectory located above the horizontal axis corresponds to the Verhulst logistic curve, widely known in biology, chemistry, sociology, and other sciences. The segment located below the horizontal axis corresponds to the evolution of aftershocks (Figure 2).

**Figure 2.** Logistic curve (on left) and aftershocks curve (on right) (second branch of logistic equation) at *X*∞ = 0.2. Dimensionless time *T* = *γt* is plotted along the horizontal axis.

The choice between the logistic and aftershock branches is made when setting the Cauchy problem for Equation (6). Evolution proceeds along the aftershock branch if the initial condition satisfies the inequality *n*(0) = *n*<sup>0</sup> > *n*∞, where *n*<sup>∞</sup> = *γ*/*σ*. Thus, in the physics of aftershocks, when setting the Cauchy problem, one should set the initial conditions under the additional constraint *n*<sup>0</sup> > *n*∞. Moreover, it is reasonable to use the strong inequality *n*<sup>0</sup> *n*<sup>∞</sup> = *γ*/*σ*. Indeed, for *t* → ∞, the frequency of aftershocks asymptotically approaches from above to the background (equilibrium) value *n*∞. Experience shows that as a rule *n*<sup>0</sup> *n*<sup>∞</sup> after a strong earthquake. The analysis of Equation (6) under the condition *n*<sup>0</sup> *n*<sup>∞</sup> indicates that at the first stages of evolution, the frequency of aftershocks decreases with time in accordance with the classical Omori Formula (1). Let us take a closer look at this important circumstance, since the existence of the aftershock branch is not so widely known.

The aftershock branch is entirely located above the saturation level *n*∞ and is a monotonically decreasing function of time. When *t* → +∞ it tends asymptotically from above the saturation level (see the right panel in Figure 2). When setting the Cauchy problem in the physics of aftershocks, the initial condition should be asked the restriction *n*<sup>0</sup> *n*∞.

Let us show that the decrease in the frequency of aftershocks with time at the first stage of evolution occurs according to the Omori hyperbola (1). It is natural to call this stage of evolution the *Omori epoch*.

Let us introduce the notation

$$t\_{\infty} = \frac{1}{\gamma} \ln(1 - \frac{n\_{\infty}}{n\_0}),\tag{8}$$

and write the solution of evolution Equation (6) in the following form:

$$n(t) = n\_{\infty} \left\{ 1 - \exp\left[\gamma (t\_{\infty} - t)\right] \right\}^{-1}.\tag{9}$$

In the Omori epoch *t*<sup>∞</sup> < *t* 1/*γ* and, respectively,

$$m(t) = 1/\sigma(t - t\_{\infty}).\tag{10}$$

Formula (10) coincides with the classical Omori Formula (1) up to notation.

Observational experience indirectly testifies to the plausibility of our logistic model. It is known, for example, that over time the frequency *n* tends not to zero, as follows from the Omori law, but to some equilibrium value *n*∞. Further, some combination of the logistic and aftershock branches makes it possible to propose a scenario for the occurrence of an earthquake swarm (see details in [11]).

#### **4. Stochastic Equation**

Changing variables in a differential equation is often a powerful tool for finding solutions to it. We already know the solutions for the Omori Equation (2) and the logistic Equation (6). Nevertheless, we will still change the variable *n*(*t*) in order to linearize both of these equations. This will make it easier for us to search for a stochastic generalization of the equation for the evolution of aftershocks.

The following replacement will help us transform nonlinear Equations (2) and (6) into linear ones [3]:

$$n(t) \to \lg(t) = 1/n(t). \tag{11}$$

Omori Equation (2) takes on an extremely simple form:

$$
\dot{\mathbf{g}} = \sigma.\tag{12}
$$

Here the dot above the symbol means time differentiation. Logistic Equation (6) becomes a first-order linear differential equation:

.

$$
\dot{\mathfrak{g}} + \gamma \mathfrak{g} = \sigma.\tag{13}
$$

Generally speaking, this circumstance is quite interesting in itself, but we use it precisely to make the stochastic generalization of the evolution equation the simplest possible way. Namely, let us imagine that the deactivation coefficient experiences small fluctuations. This assumption is formalized as follows: *σ* → *σ* + *ξ*(*t*), where *ξ*(*t*) is a random function of time, and *max*|*ξ*| *σ*. As a result, we have

$$
\dot{\mathfrak{g}} + \gamma \mathfrak{g} = \sigma + \mathfrak{f}(\mathfrak{t}).\tag{14}
$$

Let us formally solve the Equation (14):

$$\log(t) = \mathcal{g}\_{\mathcal{D}} + (\mathcal{g}\_0 - \mathcal{g}\_{\mathcal{D}}) \exp(-\gamma t) + \int\_0^t \zeta(t\_1) \exp[\gamma(t\_1 - t)] dt\_1. \tag{15}$$

Here *g*<sup>∞</sup> = *σ*/*γ*, *g*<sup>0</sup> = *g*(0).

Our second assumption is that *ξ*(*t*) is the Langevin source, i.e., delta-correlated random function with zero mean

$$
\overline{\xi(t)} = 0, \qquad \overline{\xi(t)\xi(t\_1)} = N\delta(t - t\_1), \tag{16}
$$

where the line at the top means averaging.

Now Equation (14) should be considered as the Langevin stochastic equation (e.g., see [13], where the Langevin equation is studied in detail). The phenomenological parameter N is determined by the intensity of noise affecting our dynamic system.

#### **5. Nonlinear Diffusion Equation**

If we ask ourselves how to describe the evolution of aftershocks not only in time, but also in space–time, then this immediately puts us in a difficult position. On the one hand, a number of methods are known for modeling space–time distributions, but on the other hand, in our case, there is a strong limitation, which consists in the fact that when averaging over the epicentral zone, we want to obtain the Omori law in the form (2), or in the form (6). Fortunately for us, it turns out here that we can use the well-known Kolmogorov–Petrovsky–Piskunov equation (abbreviated KPP), which describes nonlinear diffusion [14]. It is convenient for us to represent it in the following form:

$$\frac{\partial n}{\partial t} = n(\gamma - \sigma n) + D \frac{\partial^2 n}{\partial x^2} \tag{17}$$

where *n*(*x*,*t*) is the spatio-temporal distribution of aftershocks, the *x* axis is directed along the earth's surface (for simplicity, we limited ourselves to a one-dimensional model), and *D* is a new phenomenological parameter (diffusion coefficient). At *D* = 0, (17) turns into the logistic equation of the evolution of aftershocks (6), and under the additional condition *γ* = 0 into Omori law (2).

It is useful to derive (17) from the integro-differential equation

$$\frac{\partial n}{\partial t} = \Phi(n) + \int\_{-\infty}^{\infty} K(x - y) \cdot n(y, t) \, dy. \tag{18}$$

Here Φ(*n*) is some kind of functional. This will make it possible to express the parameters *γ* and *D* in terms of the kernel *K*(*x* − *y*). Indeed, suppose that *K*(*x* − *y*) = *K*(*y* − *x*), i.e., the core is symmetrical. If *K* → 0 for |*x* − *y*| → ∞, then expanding *n*(*x* − *z*, *t*) into a Taylor series in powers of *x*, we obtain

$$\frac{\partial n}{\partial t} = \gamma n + \Phi(n) + D \frac{\partial^2 n}{\partial x^2} + \dots \tag{19}$$

where

$$\gamma = \int\_{-\infty}^{\infty} K(z)dz,\ \mathrm{D} = \frac{1}{2} \int\_{-\infty}^{\infty} z^2 K(z)dz,\tag{20}$$

and *<sup>z</sup>* <sup>=</sup> *<sup>x</sup>* − *<sup>y</sup>* (e.g., see [11,15,16]). Let us <sup>Φ</sup>(*n*) = −*σn*2, and confine the first two terms in the series. In this approximation we obtain Equation (17) from which after phenomenological reduction follows the Omori law in the form (2).

When constructing a phenomenological theory of aftershocks, there is one most important condition, only one, but absolutely necessary: the phenomenological coefficients, whatever meaning we put into them, we must be able to measure experimentally, since they cannot be calculated on the basis of a more fundamental theory—because we simply do not have such a theory. In Section 6.1 we show how the deactivation factor we introduced in Section 2 is estimated. In the current section, we have significantly complicated the theory and introduced the diffusion coefficient *D*. We want to briefly describe the result of observing aftershocks, which actually led us to master equation (17), and then indicate how, at least in principle, the parameter *D* can be estimated experimentally. (In this regard, it is appropriate to mention the recent interesting work [17]. It provides additional arguments in favor of the idea of the applicability of the KPP equation for modeling aftershocks).

The study of aftershocks in space and time led us to the idea of using the KPP equation as the master equation. The main step forward in the study of the space–time distribution was the discovery that, apparently, at least some of the aftershocks tend to propagate like waves with a speed much lower than the speed of seismic waves [8,18]. The rate of propagation varies widely from case to case. It is roughly a few kilometers per hour. This value is three orders of magnitude less than the velocities of elastic waves in the crust, which suggests the propagation of a nonlinear diffusion wave excited by the main shock.

Equation (17) has self-similar solutions in the form of a traveling wave *n*(*x*, *t*) = *n*(*x* ± *Ut*) [14,19]. It is this circumstance that played a role in our choice of the KPP equation as the master equation. The estimation of the wave propagation velocity can be performed by analyzing the dimensions of the coefficients of the master equation: *<sup>U</sup>* <sup>∼</sup> <sup>√</sup>*γD*. Knowing the propagation velocity *<sup>U</sup>*, and estimating the parameter *<sup>γ</sup>* according to the formula *γ* = *n*∞*σ*, we can give an oriented estimate of the diffusion coefficient *D* = *U*2/*γ*.

#### **6. Discussion**

We have outlined a phenomenological basis, united by the general idea of a differential approach, to describing and understanding the dynamics of aftershocks. Starting with an elementary nonlinear differential Equation (2), completely equivalent to the Omori law in its classical expression (1) at *σ* = *const*, we tried to use minimal modifications in order to go first to the logistic equation (6), then to the stochastic Equation (14), and, finally, to the nonlinear diffusion Equation (17). Perhaps it would be useful to note that we explicitly used the methodological principle of Descartes, the essence of which is that one should go from simple to complex, using clear and precise modifications of the theoretical description of the problem under study.

Taken together, four phenomenological master equations, united by a common idea, make it possible to comprehend a fairly wide range of properties and patterns of aftershocks found experimentally. Moreover, phenomenological theory allows certain predictions to be made that can be verified experimentally. Let us illustrate what has been said with a number of examples.

#### *6.1. Inverse Problem*

The inverse problem of the physics of earthquake source is to determine the phenomenological coefficients from the observation data of aftershocks. Omori law in the form (2) makes it possible to formulate and solve the problem of determining the deactivation coefficient *σ* from the observation data of the frequency of aftershocks *n*. The auxiliary value *g*, which we introduced in Section 3, is conveniently written in the form of *g* = (*n*<sup>0</sup> − *n*)/*nn*0. It is easy to see that

$$
\sigma = \frac{d}{dt} \langle \mathbf{g} \rangle. \tag{21}
$$

Here, the angle brackets denote the regularization of the function *g*(*t*) calculated from observation *n*(*t*). Regularization is reduced in this case to the smoothing procedure.

Our experience indicates that the deactivation coefficient *σ* undergoes complex variations over time [8,20]. However (and this seems to us extremely important) at the first stages of evolution *σ* = *const*. The time interval during which *σ* = *const*, we called the Omori epoch. In the Omori epoch, the classical Omori law is fulfilled (1), according to which the frequency of aftershocks hyperbolically decreases over time. In our experience, the duration of the Omori epoch varies from case to case from a few days to two or three months. We noticed a tendency for the duration of the Omori epoch to increase with the increase in the magnitude of the main shock.

An interesting prediction follows from the existence of the Omori epoch, which has been reliably confirmed by experience [21]. Namely, the question of dependence of the deactivation factor on the magnitude of main shock is analyzed theoretically and experimentally. A monotonic decrease in the deactivation factor with an increase in the magnitude of the main shock, *M*0, has been reliably established. Figure 3 shows the result of measurements of *σ* at different values of *M*0. To measure the deactivation factor, we used the USGS/NEIC earthquake catalog and the technique developed during the compilation of the Atlas of Aftershocks [20]. We see that, on average, *σ* decreases monotonically with the increase in *M*0. The dependence *σ*(*M*0) is approximated by the formula

$$
\sigma = A - BM\_{0\prime} \tag{22}
$$

where *A*= 0.64, *B* = 0.07 with a sufficiently high coefficient of determination *R*<sup>2</sup> = 0.82. Thus, the theoretical inequality d*σ*/d*M*<sup>0</sup> is reliably confirmed by direct measurements. We have a wonderful harmony between theory and experiment.

**Figure 3.** Dependence of the deactivation factor of the earthquake source on the magnitude of the main shock. Other explanations see in the text below.

It is quite clear that the question of how best to formulate the Omori law, in the form (1) or (2), could only be solved by observation and experience. Equation (2) turned out to be more effective, since it made it possible to introduce a simple and useful concept of deactivation of the source, to pose the inverse problem of the source, and to reveal the existence of the Omori epoch. In addition to this, we have shown that with the help of (2) one can make a meaningful statement about the deactivation coefficient and, moreover, check this statement experimentally.

Finally, the question of whether it is not better to use the one-parameter Formula (2) for modeling aftershocks than the two-parameter Hirano–Utsu formula *n* = *k*/(*c* + *t*) *<sup>p</sup>* [2] deserves discussion. We give preference to Formula (2), since the inverse problem solved on its basis indicates the existence of the Omori epoch [8,9]. The Hirano–Utsu formula is unacceptable, since it contradicts the existence of the Omori epoch at *p* = 1, and at *p* = 1 it coincides with the Omori Formula (1).

Perhaps it would be appropriate to draw a distant historical analogy here. According to the law of gravitation, the interaction potential *ϕ* ∝ 1/*r* leads to the ellipticity of the planetary orbit. The deviation from ellipticity, for example, of the orbit of Mercury, could in no way serve as a reason for choosing the interaction, say, in the form ∝ 1/*rp*. Another understanding of the deviation of orbit from strictly elliptical had to be looked for, and it was found within the framework of general relativity. Perhaps, finding themselves in a similar situation, Hirano and Utsu should not have immediately abandoned the excellent Omori law, but should have looked for other explanations for the deviation of the real flow of aftershocks from strict hyperbolicity.

#### *6.2. Triggers*

Deviations from the classical Omori law (1) are caused not only by the nonstationarity of the parameters of the geological environment in the source, which we have expressed in the form of a possible dependence of the deactivation factor on time. Deviations can occur under the influence of so-called triggers, i.e., relatively small disturbances of geophysical fields of endogenous or exogenous origin.

We point here to two endogenous triggers that we have recently discovered. Both are aroused on the main shock. One of them has the form of a round-the-world seismic echo, and the second represents free elastic oscillations of the Earth as a whole, excited by the main shock. We have described both triggers in detail in a number of papers, so we will restrict ourselves here to references [22–26].

Let us dwell on exogenous triggers in more detail. For a long time, the cosmic effects on seismicity have been widely discussed, but there is still no agreement in the geophysical community on the effectiveness of such impacts. The controversy about the influence of geomagnetic storms on the global activity of earthquakes arises especially often (see, for example, [27–30]). The question is really difficult. On the one hand, observations indicate a correlation between seismicity and geomagnetic storms and a complex of electromagnetic phenomena associated with them (see papers [31–37] and the literature cited therein). On the other hand, the mechanism of the impact of geomagnetic storms on rocks, leading to modulation of seismicity, is not entirely clear. In this regard, the idea of the magnetoplasticity of rocks [32,33] seems to us very encouraging, but a discussion of this deep idea would lead us far astray.

A wide class of exogenous triggers of anthropogenic origin is known. We will restrict ourselves here to an indication of the weekend effect discovered in [38], and the so-called Big Ben effect, or the effect of hour markers [39,40]. Both effects pose a difficult question for the researcher about the global impact of the industrial activity of mankind on the lithosphere.

#### *6.3. Triads*

Apparently, the idea of a peculiar trinity of foreshocks, main shock, and aftershocks in a sequence of tectonic earthquakes [41–43] was formed in seismology not without the influence of mathematics, in which a binary relation between elements of a set can give rise to a trichotomic relation. The trinity of foreshocks, main shock, and aftershocks was proposed to be called the classical triad [44]. The magnitude of the main shock *M*<sup>0</sup> is always greater than the maximum magnitudes of foreshocks and aftershocks. The classical triad satisfies the inequalities

$$M\_- < M\_{+, \prime} \tag{23}$$

and

$$N\_- < N\_+. \tag{24}$$

Here *M*−(*M*+) and *N*−(*N*+) are the maximum magnitude and the number of foreshocks (aftershocks), respectively.

Quite often *N*<sup>−</sup> = 0, i.e., foreshocks are absent even before rather strong earthquakes. Figure 4 illustrates this situation (the database and the plotting method will be described in detail below). With regard to aftershocks, there is a stable opinion that after a sufficiently strong earthquake, repeated tremors are always observed, i.e., *N*<sup>+</sup> = 0.

In this section, we want to present rare but extremely interesting types of anomalous triad, for which inequalities that are directly opposite to inequalities (23) and (24) hold [45]. These are the so-called mirror triads, for which *M*<sup>−</sup> > *M*+, *N*<sup>−</sup> > *N*<sup>+</sup> and symmetric triads, for which *N*<sup>−</sup> = *N*+. Moreover, *N*<sup>+</sup> = 0 in a significant part of the mirror triads.

**Mirror triads.** Extensive literature is devoted to the experimental study of the classical triads. We point here to work [21], since in the study of anomalous triads we used a database and general methods of analysis similar to those used here in the study of classical triads (see also [9]).

**Figure 4.** Generalized picture of a shortened classical triad. Zero moment of time corresponds to the moment of the main shock.

We used data on earthquakes that occurred on Earth from 1973 to 2019 and were registered in the world USGS/NEIC catalog of earthquakes (https://earthquake.usgs.gov; last accessed on 30 April 2022). There were found *N*<sup>0</sup> = 2508 main shocks with a magnitude of *M*<sup>0</sup> ≥ 6 and a hypocenter depth not exceeding 250 km. For each main shock, a circular epicentral zone was determined by the formula *lgL* = 0.43*M*<sup>0</sup> − 1.27, where the radius of the zone *L* is expressed in kilometers [46]. According to our definition, the classical triad is formed by earthquakes, which occurred in the epicentral zone in the interval of ±24 h relative to the moment of the main shock, provided that the inequalities (23), (24) are satisfied. The total number of earthquakes was distributed among the members of the triad in the following way: *N*<sup>−</sup> = 1105, *N*<sup>0</sup> = 2398, *N*<sup>+</sup> = 31865. Note that here and below *N*<sup>0</sup> is the number of main shocks. Figure 4 is based on truncated triads, in which there are no foreshocks. The graph was constructed by the method of overlapping epochs, and the main shock of the earthquake was used as a benchmark. For truncated triads, the distribution looks like this: *N*<sup>−</sup> = 0, *N*<sup>0</sup> = 2066, *N*<sup>+</sup> = 21422. The distribution at *N*<sup>−</sup> = 0: *N*<sup>−</sup> = 1105, *N*<sup>0</sup> = 332, *N*<sup>+</sup> = 10443. We see that in the presence of foreshocks, the activity of aftershocks is higher than in the absence of foreshocks and the number of truncated triads significantly exceeds the number of complete ones.

In the course of studying classical triads, the idea arose to make a selection of observational data by replacing inequalities (23), (24) with the opposite ones. As a result, it was possible to find the mirror triads. Figure 5 shows the truncated mirror triads: *N*<sup>−</sup> = 237, *N*<sup>0</sup> = 156, *N*<sup>+</sup> = 0. If *N*<sup>+</sup> = 0, then *N*<sup>−</sup> = 1375, *N*<sup>0</sup> = 104, *N*<sup>+</sup> = 755.

**Figure 5.** Generalized view of truncated mirror triads. Zero moment of time corresponds to the moment of the main shock. The red line is obtained by averaging over a sliding window of 20 min, with the step 1 min.

We see that mirror triads are relatively rare phenomenon. They appear about an order of magnitude less frequently than classical triads. To make the picture of mirror triads more visual, we will show Figure 6. It shows mirror triads in the range of main shocks magnitude 5 ≤ *M*<sup>0</sup> < 6. Here *N*<sup>−</sup> = 4189, *N*<sup>0</sup> = 2430, *N*<sup>+</sup> = 201.

**Figure 6.** Time distribution of foreshocks and aftershocks of mirror triads in the range of magnitudes of the main shocks 5 ≤ *M*<sup>0</sup> < 6.

**GTS.** So, we found that there is a rare but rather interesting subclass of tectonic earthquakes, in which the number of aftershocks in the interval of 24 h after the main shock is significantly less than the number of foreshocks in the same interval before the main shock. In many cases, there are no aftershocks at all. We asked the question: Are there earthquakes with magnitudes *M*<sup>0</sup> ≥ 6, neither before nor after which there are neither foreshocks nor aftershocks? The search result was amazing. We have discovered a wide variety of this kind of earthquake and named it *Grande terremoto solitario* (Italian), or GTS for short [47]. In Figure 7, we see that the number of GTS (2460) is approximately equal to the number of classical triads (2398).

**Figure 7.** Solitary earthquakes with magnitudes *M*<sup>0</sup> ≥ 6.

GTS arise spontaneously under very calm seismic conditions and are not accompanied by aftershocks. This suggests an analogy between the GTS and the so-called "Rogue waves" (or "Freak waves")—isolated giant waves that occasionally emerge on a relatively quiet ocean surface (e.g., see [48]). This analogy may prove to be quite profound, since the spontaneous occurrence of pulses with anomalously high amplitudes is a common property of the nonlinear evolution of dynamic systems [49].

For completeness, we also present the data for the symmetric triads in which *M*<sup>0</sup> ≥ 6 and *N*<sup>+</sup> = *N*−: *N*<sup>−</sup> = 186, *N*<sup>0</sup> = 121, *N*<sup>+</sup> = 186. It is interesting to note that, formally, GTS can be related to a variety of symmetric triads, since for them *N*<sup>+</sup> = *N*<sup>−</sup> = 0.

**Activation factor.** Figures 5 and 6 shows that foreshocks in the mirror triad appear to have a temporal distribution similar to the Omori distribution for aftershocks in the classical triad. Let us dwell on this in more detail. We represent the classical Omori law [1] in the simplest differential form .

$$
\dot{\mathbf{g}} = \sigma\_+.\tag{25}
$$

Here *g* = 1/*n*, *n*(*t*) is the frequency of aftershocks, *t* > 0, the dot above the symbol means time differentiation, *σ*+ is the so-called deactivation factor of the earthquake source, "cooling down" after the main shock (see [3,8,21]). Suppose that for foreshocks of the mirror triad, the evolution law (25) is fulfilled with the replacement of *σ*<sup>+</sup> by *σ*−. It is natural to call the *σ*<sup>−</sup> value the activation factor.

In this paper, we will limit ourselves to presenting the interesting Figure 8. It shows the generalized evolution of foreshocks and aftershocks in symmetric triads satisfying the condition 5 ≤ *M*<sup>0</sup> < 6. Here *N*<sup>−</sup> = 1050, *N*<sup>0</sup> = 742, *N*<sup>+</sup> = 1050. The top panel shows an amazing mirror image. In the bottom panel, we have shown the variations of the *σ*<sup>−</sup> and *σ*+ functions as the first step towards studying the activation and deactivation coefficients of an earthquake source in mirror triads. (For the procedure for calculating *σ*±, see [8]).

**Figure 8.** Time dependence of the properties of symmetric triads. From top to bottom: earthquake frequency, activation (blue), and deactivation (red) factors.

**Origin of mirror triads.** In conclusion of this section, we would like, with all the necessary reservations, to express a careful judgment on the question of the origin of mirror triads. Let us assume that a system of faults in a certain volume of rocks is under the influence of a slowly growing total shear stress *τ*. Threshold tension *τ*<sup>∗</sup> at which destruction occurs, i.e., the sides of the fault shift and a rupture occurs, generally speaking, is the lower, the larger the linear dimensions of the fault *l*:

$$
\pi\_\* = \mathbb{C}l^{-m}.\tag{26}
$$

Here *C* is the dimensional coefficient of proportionality, depending on the properties of rocks in the selected volume, *m* > 0. Then, the largest fault reaches the threshold first. Its destruction is manifested in the form of the main shock of an earthquake with a magnitude *M*0.

If the *C* parameter is uniformly distributed over the source volume, then foreshocks do not arise. Aftershocks appear due to the fact that after the main shock, the general external stress is partially removed, and the remaining stress is redistributed in a complex way throughout the volume. The local over-tensions arise, in such a way that smaller faults than the one that generated the main shock can be activated and give repeated tremors. This is how we can imagine the emergence of a shortened classical triad.

In some cases, the specific distribution of faults in terms of *l* and the distribution of local stresses may turn out to be such that not a single aftershock occurs. It is possible that such a situation occurs when the GTS is excited.

The appearance of the mirror triad can be understood if we assume that the parameter *C* is not uniformly distributed over the volume, or rather, that there is a strong scatter in the *C*(*l*) values. Then a situation is possible when, before the largest fault is destroyed, smaller faults are activated, and foreshocks appear. A triad of tectonic earthquakes will appear. Whether it will be classical or mirror-like depends on the distribution of faults by the value of *l*, on the dispersion of the *C*(*l*) coefficient, and on the mosaic of local stresses that arose after the main shock.

As a summary, we point out that classical (normal) triads make up approximately 85% of all triads. Anomalous triads account for 15%, with mirrored, 10% and symmetrical, 5%. In this calculation, we excluded the GTS, which seem to form a special set of earthquakes. In the class of tectonic earthquakes, we found a subclass of the so-called mirror triads. A specific property of mirror triads is that, in contrast to classical triads, in which the number of aftershocks is greater than the number of foreshocks, in mirror triads the number of aftershocks is less than the number of foreshocks in the interval 24 h before and 24 h after the main shock. In many cases, there are no aftershocks at all. In addition to this, strong solitary earthquakes were discovered, which are not preceded by foreshocks, and after which there are no aftershocks.

The mirror triads, these ghosts of the classical triads, are not only curious in themselves, but can most decisively influence our understanding of the alternative possibilities of the dynamics of lithosphere, leading to catastrophic earthquakes. In particular, we face a fundamental problem, the essence of which is to find the physical and geotectonic reasons for the apparent predominance of truncated classical triads in the seismic activity of the Earth.

Concluding the discussion, we want to make a judgment, perhaps controversial, that Omori's law, like no other, gives us the opportunity to realize the uncertainty, incompleteness, and, in a certain sense, immaturity of the physics of earthquakes. So far, it is still possible and, in reality, there is a wide range of opinions on the unresolved issues of the Earth's seismicity. So far, everyone, in accordance with their idea of the nature of cognition, can see in the Omori law either an empirical formula convenient for approximating observations, or an indication of the deep physical meaning of hyperbolicity in the frequency of aftershocks after a strong earthquake.

#### **7. Conclusions**

We will deviate from tradition and, instead of simply listing the results of the work (which is still far from complete), we will briefly present some prospects of the research that we propose to carry out in the near future. The study will be devoted to the geometrodynamics of a tectonic earthquake source. We proceed from the fact that in the study of earthquakes some geometrically visual representations and considerations are necessary, and that analytics alone is insufficient. Representing the source as the interior of the convex hull of a point manifold of aftershocks, we have outlined a program whose purpose is to reduce the mosaic of very complex, intricate realities of the evolving source to geometric

objects. An interesting object is the space–time trajectory of the shell's center of gravity. The curvature of the envelope surface of the epicenters in dynamics and other geometrically visual images can also turn out to be quite interesting. As conceived, the geometrodynamics will become a source of new ideas for the development of the phenomenological theory of earthquakes.

However, let us return to the work presented by us and try to make some summary of all the results. The general conclusion is as follows: the methodological approach based on differential equations of evolution opens up new possibilities for the analysis of experimental material. The phenomenological equations of evolution proposed by us allows for the posing of inverse problems of the source physics and makes it possible to formulate unexpected questions regarding the dynamics of earthquakes. The phenomenological theory, a sketch of which we have given here, not only enriches the system of ideas about the source, but, we hope, indicates the possibility of searching for approaches to solving problems of a fundamental nature.

We are fully aware of the fact that neither the totality of facts, even if it is represented by a set of empirical formulas, nor a logically consistent phenomenological theory, by itself, lead us to a deep penetration into the essence of earthquakes. A deep understanding would be much more facilitated by a theory based on the fundamental laws of physics, taking into account the characteristics of the geological environment. Theoretical constructions of this kind are known, but for a completely understandable reason they refer only to individual aspects of the phenomenon, and not to the phenomenon as a whole. Under these conditions, it is reasonable and natural to consistently continue the development of the phenomenological theory of earthquakes, the foundations of which were laid by Fusakichi Omori.

**Author Contributions:** Conceptualization, A.G., O.Z., A.Z. and B.K.; data curation, O.Z. and B.K.; formal analysis, A.Z., O.Z., A.G. and B.K.; investigation, O.Z., A.G. and B.K.; methodology, O.Z., A.G. and B.K.; Validation, A.Z., O.Z., A.G. and B.K.; writing—original draft, A.G.; writing—review and editing, A.Z., O.Z., A.G. and B.K.; project administration, A.Z. and A.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Ministry of Science and Higher Education of the Russian Federation as part of the state assignment of the Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences # 075-00693-22-00.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The USGS earthquakes catalog from 1973 to 2019 is accessible from https://earthquake.usgs.gov/earthquakes/search/ (last accessed on 30 April 2022).

**Acknowledgments:** We express our deep gratitude to A.L. Buchachenko, F.Z. Feygin, N.A. Kurazhkovskaya, and A.S. Potapov for their interest in this work. The authors thank the staff U.S. Geological Survey for providing the catalogues of earthquakes.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

