Statistical Physics Model of Seismic Activation Preceding a Major Earthquake

1. Introduction

An increase in the number of intermediate sized earthquakes ($\mathrm{M}>3.5$) in a seismic region preceding the occurrence of an earthquake with magnitude $\mathrm{M}>6$, referred to as seismic activation, has been observed to occur Bowman et al. (1998). For example, seismic activation was observed in a geographic region spanning ${21}^{\circ }N-{26}^{\circ }N\times {119}^{\circ }E-{123}^{\circ }E$ for a period of time between 1991 and 1999 preceding the magnitude 7.6 Chi-Chi earthquake Chen (2003). Figure 1 shows a schematic plot of the cumulative distribution of earthquakes of different magnitudes in a region undergoing seismic activation in two different time intervals of equal duration preceding occurrence of a major ($7<\mathrm{M}<8$) earthquake at time $\tau ={\tau }_0$. In this figure, $\tau$ is a real time parameter, and ${\tau }_0$ is the characteristic time of major earthquake recurrence assuming an earthquake of similar magnitude occurred in the same region at $\tau =0$ Rundle et al. (2001); Vallianatos and Sammnds (2004). Importantly, the cumulative distribution of earthquakes in a time interval of fixed width increasingly deviates away from a Gutenberg-Richter linear log-magnitude plot as the end of the time interval approaches ${\tau }_0$.

As a means of predicting the time $\tau ={\tau }_0$ at which a major earthquake preceded by seismic activation occurs, it has been hypothesized that the average seismic moment ${\langle \mathrm{M}\rangle }_{\tau }$ of earthquakes occuring in intervals of time $(\tau ,\tau +\Delta \tau )$ preceding a major earthquake obeys an inverse power of remaining time to failure law:

$${\langle \mathrm{M}\rangle }_{\tau }\propto {\displaystyle \frac{1}{{({\tau }_0-\tau )}^{{\gamma }_1}}}$$

(1)

and that as a result, the cumulative Benioff strain $\mathcal{C}\left(\tau \right)$, defined as:

$$\mathcal{C}\left(\tau \right)=\sum _{i=1}^{n\left(\tau \right)}{\mathrm{M}}_{0,i}^{1/2},$$

(2)

where ${\mathrm{M}}_{0,i}$ is the seismic moment of the $i^{th}$ earthquake in the region starting from a time $\tau =0$ preceding the major earthquake, and $n\left(\tau \right)$ is the number of earthquakes occurring in the region up to time $\tau$, satisfies Tzanis et al. (2000):

$$\mathcal{C}\left(\tau \right)=a-b{({\tau }_0-\tau )}^{{\gamma }_2},{\gamma }_2=1-{\gamma }_1/2.$$

(3)

Notably, the exponent selection of $1/2$ in equation (2) is not necessary to derive formula (3) with a different arithmetic relation between ${\gamma }_1$ and ${\gamma }_2$, but appears to have been commonly selected by previous researchers based on resulting predictions of major earthquake occurrence time when formula (3) is fit to real seismic data Vallianatos and Chatzopoulos (2018). It is also noted that when a fit to real seismic data is performed, a value ${\gamma }_2\approx 0.3$ is typical Bowman et al. (1998).

A mathematical model of seismic activation based on damage mechanics of earthquake faults has been put forth to derive equation (3) with a value ${\gamma }_2=1/3$ Ben-Zion and Lyakhovsky (2002). In this derivation, the occurrence of seismic activation earthquakes progressively decreases the average shear modulus of fault material in the seismic region where subsequent seismic activation earthquakes occur, and the result ${\gamma }_2=1/3$ is determined by a Boltzman kinetic type description of how ruptured faults of different lengths at different positional locations grow and join together Tzanis and Vallianatos (2003).

In addition to the damage mechanics model of seismic activation, an empirical model of seismic activation using statistical physics known as the Critical Point (CP) model has been put forth to derive equation (3) with a value ${\gamma }_2=1/4$ Rundle et al. (2001). In this derivation, the inverse power of remaining time to failure law:

$${\langle \mathrm{M}\rangle }_{\tau }\propto {\displaystyle \frac{1}{{({\tau }_0-\tau )}^{3/2}}}$$

(4)

is asserted based on identifying the mean rupture length $\mathcal{L}\left(\tau \right)$ of earthquakes occuring at time $\tau$ with the correlation length of a statistical physical system described by Ginzburg-Landau mean field theory with a temperature parameter depending on $\tau$, whereby:

$$\mathcal{L}\left(\tau \right)\propto {\displaystyle \frac{1}{{({\tau }_0-\tau )}^{1/2}}},$$

(5)

and relation (4) follows from the scaling relation ${\langle \mathrm{M}\rangle }_{\tau }\propto \mathcal{L}{\left(\tau \right)}^3$ which holds when the fault material shear modulus is constant in the seismic activation region Rundle et al. (2003). Importantly, previous work has not explained why it is physically reasonable to describe statistics of seismic activation with thermal equilibrium statistical physics formalism, or what predictions other than time of major earthquake occurrence are possible with such models Newman et al. (1995). Therefore, the objective of this article is to advance the detailed mathematical description of the correspondence between nonlinear differential equation modelling and statistical physics modelling of seismic activation in a way that advances testing of model predictions against real seismic measurements.

The outline of the article is as follows. Section 2 introduces a sine-Gordon equation modelling earthquake fault dynamics during seismic activation and explains how inverse scattering theory of this equation implies a relation between solutions to a nonlinear eigenvalue problem whose eigenfunctions characterize the seismic velocity model of the region and statistical physics. Section 3 explains how statistics of the eigenvalue solutions are related to earthquake occurrence statistics, and how describing these statistics with a 2D Coulomb gas statistical physics model accounts for deviation of earthquake occurrence statistics from Gutenberg-Richter statistics before a major earthquake. Section 4 concludes by commenting on how statistical physical modelling implies the phase space dimension of a nonlinear dynamical system characterizing rupture nucleation in the seismic activation region is finite, and how this implication can be tested against real seismic measurements.

2. Methods

2.1. 1D Fault Dynamics Inverse Scattering Theory

In 1+1 spacetime dimensions, the differential equation:

$$A{\partial }_{\tau }^{2}U(\tau ,z)+B{\partial }_{\tau }U(\tau ,z)-C{\partial }_z^{2}U(\tau ,z)=-sin(U(\tau ,z)/D).$$

(6)

has been used to model migration of earthquake hypocentres along earthquake faults in seismic regions over periods of time during which multiple earthquakes occur Bykov (2001). In this equation, $\tau$ is real time, z coordinates a direction of earthquake hypocenter migration along an earthquake fault, $U(\tau ,z)$ is the local displacement of elastic material across the earthquake fault, $A{\partial }_{\tau }^{2}U(\tau ,z)$ is the local inertial force acting on the fault material, $B{\partial }_z^{2}U(\tau ,z)$ is the local elastic restoring force acting on the fault material, and $C{\partial }_{\tau }U(\tau ,z)$ and $sin\left(U\right(\tau ,z)/D)$ are local frictional force acting on the fault material attributed to periodic contact of the material with tectonic plates on either side of the fault. If the earthquake fault material has constant height h and shear modulus $\mu$ along the fault, a soliton (2 kink) solution to equation (6) can be interpreted to describe growth of an interval where earthquake rupture has occurred along the z-axis, which in this 1 dimensional model is identified with the both the direction of hypocentre migration and the direction of all earthquake slip vectors Carlson et al. (1994).

Restricting focus to the case $C=0$, with rescaling of $\tau$, z, and $U(\tau ,z)$, each of the constants A, B, and D in equation (5) can be set to 1. With this rescaling, and definition of the matrices:

$$\begin{array}{ccc}\hfill \mathcal{M}& =& \left[\begin{array}{cc}-i\omega & -{\displaystyle \frac{1}{2}}{U}_z(\tau ,z)\\ {\displaystyle \frac{1}{2}}{U}_z(\tau ,z)& i\omega \end{array}\right],\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \mathcal{N}& =& {\displaystyle \frac{i}{4\omega }}\left[\begin{array}{cc}cosU(\tau ,z)& sinU(\tau ,z)\\ sinU(\tau ,z)& -cosU(\tau ,z)\end{array}\right],\hfill \end{array}$$

(8)

for an arbitary complex number $\omega$, the equation:

$${\mathcal{M}}_{\tau }-{\mathcal{N}}_z+\mathcal{M}\mathcal{N}-\mathcal{N}\mathcal{M}=0,$$

(9)

is equivalent to equation (5) Khan (2020). Introduction of the matrix operator $\mathcal{M}$ also permits introduction of the associated seismic wave scattering problem:

$${\displaystyle \frac{\partial }{\partial z}}\left[\begin{array}{c}{\Psi }_1(\tau ,z)\\ {\Psi }_2(\tau ,z)\end{array}\right]=\mathcal{M}\left[\begin{array}{c}{\Psi }_1(\tau ,z)\\ {\Psi }_2(\tau ,z)\end{array}\right],$$

(10)

in which $\tau$ is an auxillary parameter, and oscillatory dependence $e^{-i\omega t}$ of the incoming and outgoing scattered waves on a scattering time parameter t has been implicitly assumed in the definitions of $\mathcal{M}$ and $\mathcal{N}$. According to the inverse scattering method, the importance of the operator $\mathcal{N}$ is that it determines real time $\tau$ evolution of scattering data determined by the operator $\mathcal{M}$.

Upon definition of scattering problem (10), the inverse scattering method proceeds by determining an infinite set of its left and right scattering (i.e. Jost) solutions ${\Psi }_{\lambda ,L}(\tau ,z)$ and ${\Psi }_{\lambda ,R}(\tau ,z)$, indexed by complex wave numbers $\lambda =\omega$, with asymptotics Aktosun et al. (2010):

$$\begin{array}{ccc}\hfill {\Psi }_{\lambda ,L}(\tau ,z)& =& \left[\begin{array}{c}0\\ e^{i\lambda z}\end{array}\right],z\to \infty \hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\Psi }_{\lambda ,L}(\tau ,z)& =& \left[\begin{array}{c}{\displaystyle \frac{L(\lambda ,\tau )e^{-i\lambda z}}{T\left(\lambda \right)}}\\ {\displaystyle \frac{e^{i\lambda z}}{T\left(\lambda \right)}}\end{array}\right],z\to -\infty ,\hfill \end{array}$$

(12)

and:

$$\begin{array}{ccc}\hfill {\Psi }_{\lambda ,R}(\tau ,z)& =& \left[\begin{array}{c}{e}^{-i\lambda z}\\ 0\end{array}\right],z\to -\infty \hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill {\Psi }_{\lambda ,R}(\tau ,z)& =& \left[\begin{array}{c}{\displaystyle \frac{e^{-i\lambda z}}{T\left(\lambda \right)}}\\ {\displaystyle \frac{R(\lambda ,\tau )e^{i\lambda z}}{T\left(\lambda \right)}}\end{array}\right],z\to \infty .\hfill \end{array}$$

(14)

Next, using equation (9), real time $\tau$ evolution of the Jost functions is calculated in terms of $\tau$-variation of the reflection coefficients $R(\lambda ,\tau )$ and a finite set of complex numbers $\left\{c_j\right\}$ associated with complex wave numbers $\left\{{\lambda }_j\right\}$ of bound state solutions to linear system (10). Finally, the real time evolved sine-Gordon equation solution $U(\tau ,z)$ is determined by applying an inverse scattering transform to the real time evolved scattering data as indicated in Figure 2. Note that according to 1D scattering theory, scattering bound states are in correspondence with zeroes of the function $T\left(\lambda \right)$, while resonant scattering states are in correspondence with zeros of the reflection coefficients $R(\lambda ,\tau )$ and $L(\lambda ,\tau )$ at fixed values of $\tau$. Also note that in general the bound and resonant state $\left\{{\lambda }_j\right\}$ values have non-zero imaginary components, and are located symmetrically with respect to the imaginary axis in the complex $\lambda$-plane.

To clarify the physical relation between scattering data used by the inverse scattering method and seismic activation physics, consider the simpler 1D scattering problem in which a scattering potential function $V\left(z\right)$ is compactly supported along the z-axis, and the operator:

$$-B{\partial }_z^2+V\left(z\right),$$

(15)

has infinitely many resonant state eigenfunctions $\Psi \left(z\right)$ satisfying the elastic wave equation:

$$-B{\partial }_z^2\Psi \left(z\right)+V\left(z\right)\Psi \left(z\right)=E\Psi \left(z\right),$$

(16)

and finitely many bound state eigenfunctions ${\Psi }_j\left(z\right)$ with negative real eigenvalues $E_j={\omega }_j^2<0$. In this situation, the Jost eigenfunctions of operator (16) define a spectral kernel whose determinant equates to the potential V, which in a 1D scattering context implies the Jost functions define the velocity model of the 1D elastic wave equation Dyson (1976). Based on this established mathematical result, it is conjectured that the scattering potential function $U_z(\tau ,z)$ also equates to a determinant of an operator defined by Jost functions of scattering problem (10).

To clarify the physical relation between inverse scattering theory bound states and seismic activation it is noted that a solution $\overline{\Psi }(t,\tau ,z)$ to the linear seismic wave equation:

$${\partial }_t^2\overline{\Psi }(t,\tau ,z)-B\left(\tau \right){\partial }_z^2\overline{\Psi }(t,\tau ,z)+V(\tau ,z)\overline{\Psi }(t,\tau ,z)=0,$$

(17)

in which the auxillary seismic wave scattering time parameter t is introduced and $\tau$ tracks scattering potential changes associated with seismic activation, and for which $\overline{\Psi }(0,\tau ,z)$ and ${\partial }_t\overline{\Psi }{(t,\tau ,z)|}_{t=0}$ are compactly supported, has a resonant scattering expansion of the form Bindel and Zworski (2007):

$$\begin{array}{ccc}\hfill \overline{\Psi }(t,\tau ,z)=\sum _{j=1}^N{e}^{t\sqrt{-E_j}}{a}_j{\Psi }_{j,a}(\tau ,z)& +& e^{-t\sqrt{-E_j}}{b}_j{\Psi }_{j,b}(\tau ,z)\hfill \\ & +& \sum _{\omega resonant}{e}^{-it\omega }{c}_{\omega }{\Psi }_{\omega }(\tau ,z),\hfill \end{array}$$

(18)

in the limit $t\to \infty$. In this expansion, bound state eigenfunctions determine characteristic length scales at which unstable growth of fault material displacement may occur across the earthquake fault. Based on this interpretation, it is hereby suggested that inverse scattering thery bound states describe patterns of stress distribution in unruptured material as fault ruptures grow and nucleate.

For example, if $V(\tau ,z)$ is a potential well of depth $V_0>0$ which is independent of $\tau$, nonzero for $\left|z\right|\le L$, and zero elsewhere, there exist finitely many bound state eigenfunctions ${\Psi }_j\left(z\right)$ which decay exponentially with increasing $\left|z\right|$:

$${\Psi }_j\left(z\right)\propto e^{-k_j\left|z\right|},\left|z\right|\to \infty$$

(19)

for a discrete set of spatial decay constants:

$$k_j=\sqrt{-E_j/B},$$

(20)

whose inverse values determine characteristic length scales at which unstable growth of fault material displacement may occur across the earthquake fault. Figure 3 shows a plot of resonant and bound state frequency locations for three situations in which $-{\displaystyle \frac{1}{B}}V\left(z\right)$ is a square well potential of increasing width and fixed height.

2.2. 3D Fault Dynamics to Statistical Physics

To generalize the previous discussion of inverse scattering theory in 1 spatial dimension to 3 spatial dimensions, suppose a major earthquake hypocentre resides in an elastic half space in such a way that the elastic parameters of the half space are constant outside a sphere of diamater ${\mathcal{L}}_0$ centred at the hypocentre. Moreover, suppose that with appropriate definition of a perfectly matched layer at the sphere’s boundary, dependent on complex frequency, elastic scattering resonant frequencies at time $\tau$ are determined by solution to a nonlinear eigenvalue problem whose solutions satisfy an equation Bindel (2006):

$$R(\omega ,\tau )=0.$$

(21)

With this supposition, the complex resonant frequency set $\left\{{\omega }_j\right\}$ constitutes a 3D generalization of inverse scattering data that depends on $\tau$, and whose real time evolution is necessarily in correspondence with real time evolution of the half space elastic parameters.

Importantly, the discussion in the previous paragraph suggests that rather than being a set of continuously changing values, the resonant frequency set $\left\{{\omega }_j\right\}$ is discontinuously updated with occurrence of each seismic activation earthquake, whereby a frequency is removed from the set and other frequencies in the set are affected by alteration of the fault rupture pattern, as illustrated in Figure 4 Lei and Sornette (2022). For this reason, rather than describing the resonant frequency set deterministically, it is now put forth that this set be described probabilistically by defining a resonant frequency probability density that depends on time. More specifically, it is suggested that the resonant frequency condition $R(\omega ,\tau )=0$ be interpreted as an eigenvalue condition on a matrix selected from a $\tau$-dependent ensemble of random matrices. From this point of view, varying the real time parameter $\tau$ varies the statistical distribution of the random matrices, and in accordance with original statistical physics models of seismic activation, these random matrix models may be related by renormalization group flow Saleur et al. (1996).

3. Results

To relate the discussion in the previous chapter to seismic activation earthquake occurrence statistics, first define ${\omega }_c\left(\tau \right)$ and $\mathcal{L}\left(\tau \right)$ as the corner frequency and rupture length of the largest seismic activation earthquake occurring before time $\tau <{\tau }_0$. From this definition, the observed scaling relation between earthquake corner frequency and seismic moment $M_s$:

$${\omega }_c\left(\tau \right)\propto M_s^{-1/3},$$

(22)

together with the scaling relation:

$$M_s\propto \mathcal{L}{\left(\tau \right)}^3,$$

(23)

imply Aki (1967):

$$\mathcal{L}\left(\tau \right)\propto {\omega }_c{\left(\tau \right)}^{-1}.$$

(24)

Next, suppose that the eigenfunction of a solution $\omega$ to $R(\omega ,\tau )=0$ defines a lithospheric stress distribution preceding a seismic activation earthquake of rupture length $\mathcal{L}\propto {\omega }^{-1}$, and that the probability of this earthquake occurring the time interval $(\tau ,\tau +\Delta \tau )$ is proportional to $\omega d\tau /Q$, where $Q\left(\omega \right)$ is the quality factor of the lithospheric resonance at frequency $\omega$. Then, if $\rho \left(\omega \right)$ is the density of nonlinear eigenvalue solutions in the interval $(\omega ,\omega +d\omega )$, the number of earthquakes with corner frequency less than or equal to $\overline{\omega }$ occurring during the time interval is:

$${\dot{N}}_{c}d\tau ={\int }_{{\omega }_c\left(\tau \right)}^{\overline{\omega }}(\omega d\tau /Q\left(\omega \right))\rho \left(\omega \right)d\omega$$

(25)

To specify the mathematical form of the integral in equation (25), recall that the Gutenberg-Richter law implies the total number of earthquakes of Richter magnitude in the interval $(M_R,M_R+d{M}_R)$ occurring in the seismic activation region in the time interval $(\tau ,\tau +d\tau )$ is proportional to:

$${10}^{-b{M}_R}d{M}_{R}d\tau ,$$

(26)

which according to the relation:

$$M_R=\left({log}_{10}\left(M_s\right)-9\right)/1.5,$$

(27)

between Richter magnitude and seismic moment satisfies:

$${10}^{-b{M}_R}d{M}_{R}d\tau \propto M_s^{-1-b/1.5}d{M}_{s}d\tau \propto {\omega }^{2b-1}d\omega d\tau .$$

(28)

Therefore, assuming the Gutenberg-Richter law is valid, it follows that:

$$\rho \left(\omega \right)/Q\left(\omega \right)\propto {\omega }^{2b-2}.$$

(29)

To account for modification to the Gutenberg-Richter law in time intervals preceding a major earthquake using statistical physics, it is now conjectured that for earthquake corner frequencies:

$$\omega \approx {\omega }_c\left({\tau }_0\right)\equiv {\omega }_0,$$

(30)

the quantity $Q\left(\omega \right)/\rho \left(\omega \right)$ describes the density of a 2D Coulomb gas in the vicinity of a test charge located at $\omega ={\omega }_0$, whereby Goldenfeld (1992):

$$Q\left(\omega \right)/\rho \left(\omega \right)\propto {(\omega -{\omega }_0)}^{{\beta }_0},{\beta }_0>0,$$

(31)

for ${\omega }_0$ equal to the lowest corner frequency of a seismic activation earthquake preceding the major earthquake. With this conjecture, equation (25), modified to account for occurrence of an earthquake at corner frequency ${\omega }_0$, implies:

$${\dot{N}}_c\propto 1+{\int }_{{\omega }_0}^{\overline{\omega }}\omega {(\omega -{\omega }_0)}^{-{\beta }_0}d\omega ,$$

(32)

which in turn implies the logarithm:

$${log}_{10}{\dot{N}}_c\approx {log}_{10}\left(1+c{(\omega -{\omega }_0)}^{1-{\beta }_0}\right),$$

(33)

when plotted against Richter magnitude $M_R\propto -2{log}_{10}\omega$ for ${\beta }_0<1$, can have either of the cumulative distribution curve shapes shown in Figure 1 for different time intervals, depending on the value of ${\beta }_0$.

Having conjectured 2D Coulomb gas statistical physics are relevant to accounting for deviation of earthquake occurrence statistics from Gutenberg-Richter statistics during seismic activation, it is now further conjectured, in accordance with previous statistical physics models of seismic activation, that ${\beta }_0=\beta \left({\tau }_0\right)$, where $\beta \left(\tau \right)$ is a parameter in a $\tau$-dependent 2D sine Gordon statistical physics model whose parameters at different values of $\tau$ are related by renormalization group flow Balog et al. (2018); Carpentier (1999). With this additional conjecture, an increase in the value of $\beta \left(\tau \right)$ as $\tau \to {\tau }_0$ accounts for increasing steepness of the cumulative distribution curve shown in Figure 1. A phase diagram of a 2D sine Gordon model with a renormalization group flow indicated by arrows is shown in Figure 5. In this diagram, different vaues of the flow coordinate `t’ equate to $\beta \left(\tau \right)$ at different values of $\tau$ in such a way that $\beta \left({\tau }_0\right)$ is the horizontal coordinate of a point of tangency between two of the Ford circles. It is noted that this phase diagram corresponds to the $k=1$ description of a more general 2D sine Gordon statistical physics model defined by a field theory with $4k$ fields Dubrovin and Yang; (2020) Varchenko (1990); Zabrodin (2010).

4. Discussion

Previous work has identified predicting the time of occurrence of major earthquakes as a possible application of statistical physics models of seismic activation, but this application has not yet been realized Bowman et al. (1998). In more recent times, earthquake early warning algorithms such as FinDer and Virtual Seismologist have been developed which can in principle use previous earthquake occurrence statistics as input, and most recently, artificial intelligence algorithms such as QuakeGPT have been developed for predicting the occurrence of major earthquakes using seismic event records created with stochastic simulators as training data Böse et al. (2023); Rundle et al. (2024). Therefore, a practical applied science goal for the statistical physics model presented in this article appears to be improving the predictive performance of one or more of these existing earthquake early warning algorithms by appropriately modifying their earthquake occurrence statistical inputs, acknowledging that preliminary tests of the model’s validity against real seismic data must be passed before achieving this application objective can be considered a realistic possibility.

From a geophysical testing point of view, if it is true that the growth of unstable stress modes within the Earth during seismic activation are determined by statistical physics renormalization group flow mathematics, and, as a result, a nonlinear dynamical system of phase space dimension N characterizes the nucleation of shear stress in a seismic region preceding a major earthquake, a geophysical signal processing technique known as singular spectrum analysis should apply to determine this phase space dimension Broomhead and King (1986). Therefore, it is suggested that coda wave interferometry measurements of relative changes in seismic surface wave and/or body wave velocity be performed between pairs of seismic stations in a seismic region over a duration of time during which seismic activation is known to have occurred, and used as input to a time domain multichannel singular spectrum analysis algorithm Merrill et al. (2023). The number of channels of this algorithm would equate to the number of station pairs, and the number of singular values output by the algorithm in different time windows preceding occurrence of a major earthquake should provide some indication of a finite value for N if the statistical physics model of seismic activation is correct in principle. With reference to previous geophysical application of singular spectrum analysis, performed in the frequency domain, the signal processing algorithm suggested here is different in that it should be carried out in the time domain $\tau$ rather than the frequency domain Sacchi (2009).

In conclusion, work towards improving current earthquake early warning systems can proceed in two directions. Firstly, as an initial check on whether or not the statistical physics modelling approach presented here could be of practical utility, work can be done to determine whether or not changes of the Earth’s elastic velocity model preceding major earthquakes, as determined by coda wave interferometry, can be processed to extract an integer identifiable as the phase space dimension of a nonlinear dynamical system. Secondly, work can be done to elaborate upon the statistical physics mathematical model of seismic activation presented in this article to determine other tests of its scientific validity and potential for practical application.