Seismic Activation Modeling with Statistical Physics

1. Introduction

An increase in the number of intermediate sized earthquakes of moment magnitude $M>3.5$ in a seismic region preceding the occurrence of a mainshock event, referred to as seismic activation, has been documented by various researchers [5]. 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 [9]. Figure 1 shows a schematic plot of the cumulative distribution of earthquakes of different magnitudes in a seismic activation region in two different time intervals of equal duration preceding occurrence of a major ($7<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$[21,30]. Importantly, the cumulative distribution of earthquake magnitudes in a time interval of fixed width increasingly deviates away from a Gutenberg-Richter distribution as the end of the time interval approaches ${\tau }_0$.

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

$${\overline{\mathcal{M}}}_{\tau }\propto {\displaystyle \frac{1}{{({\tau }_0-\tau )}^{{\gamma }_1}}}$$

(1)

and that the cumulative Benioff strain $C\left(\tau \right)$, defined as:

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

(2)

where ${\mathcal{M}}_{0,i}$ is the seismic moment of the $i^{th}$ earthquake in the region starting from a time $\tau =0$ preceding the mainshock event, and $n\left(\tau \right)$ is the number of earthquakes occurring in the region up to time $\tau$, satisfies [28]:

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

(3)

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 selected by previous researchers based on Benioff’s finding that the elastic rebound of an earthquake is proportional to the square root of its seismic moment. When formula (3) is fit to real seismic data, a typical value of ${\gamma }_2$ is $0.3$ [5,29]. Notably, validity of Equation (1) has been questioned by some researchers who claim measurements of seismic activation can be explained in terms of mainshock event foreshock occurrence without acceleration of seismic release [14,32].

A model of seismic activation based on fault damage mechanics (FDM) has been used to derive Equation (3) with a value ${\gamma }_2=1/3$[3]. In this derivation, the seismic activation region is modeled as a 1D fault consisting of perfectly plastic material whose effective shear modulus decreases to 0 as $\tau \to {\tau }_0$. The result ${\gamma }_2=1/3$ is obtained from an effective shear modulus evolution equation based on non-equilibrium thermodynamic considerations [2].

In addition to the FDM model of seismic activation, a statistical mechanics model of seismic activation known as the Critical Point (CP) model has been put forth to derive Equation (3) with a value ${\gamma }_2=1/4$ [21]. In this derivation, the inverse power of remaining time to failure law:

$${\overline{\mathcal{M}}}_{\tau }\propto {\displaystyle \frac{1}{{({\tau }_0-\tau )}^{3/2}}}$$

(4)

is asserted based on identifying the mean rupture length ${\overline{\mathcal{R}}}_{\tau }$ of earthquakes occuring at time $\tau$ with the correlation length of a statistical mechanical system described by a time dependent Ginzburg-Landau equation, whereby:

$${\overline{\mathcal{R}}}_{\tau }\propto {\displaystyle \frac{1}{{({\tau }_0-\tau )}^{1/2}}},$$

(5)

which implies relation (4) assuming ${\overline{\mathcal{M}}}_{\tau }\propto {\overline{\mathcal{R}}}_{\tau }^3$ [22]. Table 1 shows typical fault material displacements and rupture lengths for earthquakes of different moment magnitudes.

Importantly, previous work on the CP model has not explained why it is physically reasonable to describe seismic activation earthquake occurrence statistics with thermal equilibrium statistical mechanics formalism. Therefore, the first objective of this article is to conjecture how FDM and CP models of seismic activation are in correspondence with each other and provide computational evidence for this conjecture. The second objective is to outline how the conjectured correspondence can be tested against seismic measurements.

To meet these objectives, two different pre-existing seismicity models are developed: the earthquake cascade model and the renormalization group theory of earthquakes [25,27]. The earthquake cascade model provides an explanation for earthquake occurrence statistics in windows of time preceding a mainshock event in terms of the statistics of metastable clusters of fault material of different size within the activation region that can join together to create larger clusters or undergo slip events to initiate earthquakes. The renormalization group theory of earthquakes explains how the formalism of equilibrium statistical mechanics accounts for time scaling of cumulative Benioff strain during seismic activation in terms of statistical mechanics critical scaling theory, without explicitly identifying the relevant statistical mechanical system. For the purposes of this article, the earthquake cascade model is developed into a random matrix model of the seismic activation region stiffness matrix, and the renormalization group theory of earthquakes is developed into a description of the time evolution of the random stiffness matrix model. In so doing, correspondence between FDM and CP seismic activation models is clarified, because the random stiffness matrix is an FDM description of fault material elastoplastic deformation, and time evolution of the random stiffness matrix is specified by statistical mechanics critical scaling theory. The conjectured correspondence between FDM and CP seismic activation models is of applied science interest because it implies a finite dimensional nonlinear dynamical system governs real time evolution of the activation region elastic model before mainshock event occurrence, and this implication can be tested against seismic measurements.

The outline of the article is as follows. Section 2 introduces the earthquake cascade model with reference to laboratory studies of rock fracture, and explains how this model is related to a random matrix model of the seismic activation region stiffness matrix. Section 3 provides a computation of metastable cluster size distribution in terms of random matrix eigenvalue statistics, conjectures how time evolution of the random stiffness matrix model is described by 2D statistical mechanics critical scaling theory, and outlines how this conjecture can be tested against seismic measurements. Section 4 comments on application of seismic activation statistical physics models.

2. Materials and Methods

2.1. Earthquake Cascade Model

It has been observed that fracture processes occurring within the Earth preceding unstable slip along a mainshock fault are analogous to fracture processes occurring in triaxially loaded rock samples preceding sample failure [19]. Triaxial test results indicate 4 stages of rock deformation leading to sample failure [10]:

• crack closure
• elastic deformation
• stable microcrack growth and extension
• unstable crack extension

Triaxial loading of sandstone samples with recording of acoustic emissions has further indicated that before unstable crack extension occurs, microcracks nucleate into clusters along faults where unstable cracks later propagate, in a process called shear localization [11]. These clusters are metastable in that they gradually grow in size before unstable crack extension occurs rapidly, and are the laboratory analog of clusters of metastable earthquake fault material whose size distribution is quantified by the earthquake cascade model.

The earthquake cascade model describes earthquake occurrence along a 2D earthquake fault. It defines $N_n\left(\tau \right)$ as the number of metastable clusters of fault material of area $2^n$ at time $\tau$, where a metastable cluster is defined as the region over which an earthquake slip propagates after initiation, and provides time evolution equations for these numbers based on cluster joining and earthquake slip initiation rates. Depending on selection of these rates, different earthquake occurrence statistics are output by the model. For example, at sufficiently low rates of earthquake occurrence, time invariant metastable cluster size statistics are:

$$N_2={\displaystyle \frac{N_1}{2^{b^{\prime }}}},N_3={\displaystyle \frac{N_2}{2^{b^{\prime }}}},...,N_n={\displaystyle \frac{N_{n-1}}{2^{b^{\prime }}}},$$

(6)

for some constant $b^{\prime }$, whereby:

$$N_k=N_1{\left(2^{-b^{\prime }}\right)}^{k-1},$$

(7)

and:

$$2^k{N}_k\propto {\left(2^k\right)}^{1-b^{\prime }}.$$

(8)

Under the assumption earthquake occurrence rate is proportional to metastable cluster area, Equation (8) specifies a Gutenberg-Richter scaling relation between frequency of earthquake occurrence and earthquake rupture area with b-value $b=b^{\prime }-0.5$. To explain this statement, note that if $N_{cM}$ is the number of earthquakes of moment magnitude greater than or equal to M occurring during the time interval $(\tau ,\tau +\Delta \tau )$, the Gutenberg-Richter law implies:

$$N_{cM}\propto {10}^{-bM},$$

(9)

which using the seismic moment-magnitude relation:

$$M=\left({\displaystyle \frac{{log}_{10}\mathcal{M}-9}{1.5}}\right),$$

(10)

and corner frequency scaling relation:

$$\mathcal{M}\propto {\omega }^{-3},$$

(11)

implies the frequency $f\left(\omega \right)$ of earthquakes with corner frequency $\omega$ and rupture area proportional to ${\omega }^{-2}$ satisfies:

$$f\left(\omega \right)\propto {\omega }^{2b-1}={\left({\omega }^{-2}\right)}^{0.5-b}.$$

(12)

To develop the earthquake cascade model to account for a continuum of possible metastable cluster sizes, the approach presented here is to modify a stochastic process model of landslide occurrence to model earthquake occurrence. Importantly, the details of the landslide model depend on whether the landslide is in a condition of secondary or tertiary creep, so its modification to describing seismicity depends on whether cumulative Benioff strain is increasing at a constant rate or accelerating.

2.2. FDM Random Matrix Model

A scalar Kesten process is a stochastic process:

$$v_{k+1}={\lambda }_{k+1}{v}_k+{ϵ}_{k+1},$$

(13)

where k is a discrete time index, $v_0=1$, and ${\lambda }_k$ and ${ϵ}_k$ are identical and independently distributed random variables. This process has been used to model sequences of landslide velocity measurements during secondary creep by assuming ${\lambda }_k$ and ${ϵ}_k$ are log-normally distributed with $\langle ln{\lambda }_k\rangle <0$, in which case the sequence of velocities are distributed according to an inverse gamma distribution $P\left(v\right)dv$ satisfying:

$$P\left(v\right)\propto v^{-b^{\prime \prime }-1},$$

(14)

for some constant $b^{\prime \prime }$ for sufficiently large velocities v [16].

In modifying process (13) to model earthquake occurrence during periods of non-accelerating cumulative Benioff strain, it is logical that the seismic stochastic process reflect classical mechanical time evolution of the seismic region, and output a random variable whose distribution accounts for the Gutenberg-Richter law. To this end, rather than describing time evolution of a single landslide velocity, it is suggested that the seismic stochastic process describe the time evolution of velocities of slow strain waves passing through the seismic activation region, and thereby also describe time evolution of metastable clusters of fault material formed by coherent motion of two or more slow strain waves. In 1 spatial dimension, slow strain waves propagating along crustal faults have been analyzed as kink solutions to a sine-Gordon equation with solitary wave velocity profiles:

$$v(x,\tau )\propto Vsech\left({\displaystyle \frac{x-V\tau }{\delta }}\right),$$

(15)

where the wave velocity V can range from a velocity of seismicity migration (${10}^{-3}$ to ${10}^{-2}$ m/s) to an earthquake rupture velocity (${10}^2$ to ${10}^3$ m/s), and $\delta$ is the thickness of fault material [7]. In terms of inverse scattering theory of the sine-Gordon equation, the velocity and width of kink solutions satisfy:

$$\begin{array}{ccc}\hfill V& =& \hfill V_r\left({\displaystyle \frac{4{\kappa }^2-1}{4{\kappa }^2+1}}\right),\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill \delta & =& \hfill {\displaystyle \frac{{\delta }_0}{\kappa +\frac{1}{4\kappa }}},\hfill \end{array}$$

(17)

where $V_r$ is the earthquake rupture velocity, ${\delta }_0$ is a characteristic fault material thickness, and $\kappa =-i\omega {\delta }_0/V_r$ is a complex number that is real when the frequency $\omega$ is a pure imaginary number identifying a 1D bound state of the sine-Gordon scattering problem [1].

The inverse scattering description of solitary waves in 1 spatial dimension can be generalized to 3 spatial dimensions using finite element analysis. Specifically, suppose the seismic activation region is contained within a hemisphere $\mathcal{H}$ centered on the Earth’s surface, and subsurface material within the hemisphere is meshed with finitely many structural elements. A nodal displacement vector $\mathbf{U}\left(t\right)=e^{-i\omega t}\mathbf{u}$ perturbing a finite strain elastoplastic deformation of material within the hemisphere at time $\tau +t$ then satisfies an equation:

$$\left(\tilde{K}(\omega ;\tau )-i\omega \tilde{D}(\omega ;\tau )-{\omega }^2\tilde{M}\right)\mathbf{u}=\mathbf{F}(\omega ;\tau )$$

(18)

where $\tilde{M}$, $\tilde{D}(\omega ;\tau )$, and $\tilde{K}(\omega ;\tau )$ are the frequency dependent finite element mass, damping, and tangent stiffness matrices at time $\tau$, and $\mathbf{F}(\omega ;\tau )$ is the vector of external forces acting on the nodes [4]. Assuming the external force is 0, the condition:

$$det\left(\tilde{K}(\omega ;\tau )-i\omega \tilde{D}(\omega ;\tau )-{\omega }^2\tilde{M}\right)=0,$$

(19)

may specify finitely many pure imaginary eigenvalues $\omega$ and associated eigenvectors $\mathbf{u}$ characterizing bound states in 3 spatial dimensions.

When the stiffness and damping matrices and external forcing vector can be approximated as being independent of $\omega$, Equation (18) can be written:

$$-i\omega \left[\begin{array}{c}\mathbf{u}\\ -i\omega \tilde{M}\mathbf{u}\end{array}\right]=\left[\begin{array}{cc}0& {\tilde{M}}^{-1}\\ -\tilde{K}\left(\tau \right)& -\tilde{D}\left(\tau \right){\tilde{M}}^{-1}\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ -i\omega \tilde{M}\mathbf{u}\end{array}\right]+\left[\begin{array}{c}0\\ \mathbf{F}\left(\tau \right)\end{array}\right],$$

(20)

which in absence of external forcing, Equation (20) implies $-i\omega$ is an eigenvalue of the matrix:

$$\tilde{L}\left(\tau \right)=\left[\begin{array}{cc}0& {\tilde{M}}^{-1}\\ -\tilde{K}\left(\tau \right)& -\tilde{D}\left(\tau \right){\tilde{M}}^{-1}\end{array}\right].$$

(21)

Therefore, if there exists an operator $\tilde{P}\left(\tau \right)$ forming a Lax pair with $\tilde{L}\left(\tau \right)$:

$${\displaystyle \frac{d\tilde{L}\left(\tau \right)}{d\tau }}=[\tilde{P}\left(\tau \right),\tilde{L}\left(\tau \right)]$$

(22)

the bound state spectrum (i.e. real eigenvalues) of $\tilde{L}\left(\tau \right)$ remains constant as $\tau$ increases, while the eigenvectors of $\tilde{L}\left(\tau \right)$ depend on $\tau$, with the equation:

$${\displaystyle \frac{d}{d\tau }}\left[\begin{array}{c}\mathbf{u}\left(\tau \right)\\ -i\omega \tilde{M}\mathbf{u}\left(\tau \right)\end{array}\right]=P\left(\tau \right)\left[\begin{array}{c}\mathbf{u}\left(\tau \right)\\ -i\omega \tilde{M}\mathbf{u}\left(\tau \right)\end{array}\right]$$

(23)

specifying this time dependence.

To pass from a deterministic to stochastic description of seismic activation region elastodynamics, differential Equation (22) may be modified to a discrete time stochastic process that outputs a sample path for $\tilde{L}\left(\tau \right)$, now regarded as belonging to an ensemble of random matrices. From this point of view, the bound state eigenvalues of the operator $\tilde{L}\left(\tau \right)$, indexed ${\omega }_1,{\omega }_2,...,{\omega }_m$, may be distributed according to a random matrix eigenvalue distribution:

$$P({\omega }_1,{\omega }_2,...,{\omega }_m)=\prod _{k=1}^{m}F({\omega }_k,\beta ){\left|\Delta ({\omega }_1,{\omega }_2,...,{\omega }_m)\right|}^{\beta },$$

(24)

where $\Delta ({\omega }_1,{\omega }_2,...,{\omega }_m)$ is the Vandermonde determinant of the bound state eigenvalues, $F(\omega ,\beta )$ is a confining potential of the bound state eigenvalues in 1 spatial dimension, and $\beta$ is henceforth referred to as a damage parameter [12]. Physically, in analogy to dependence of kink velocity on bound state eigenvalue expressed by Equation (16), distribution (24) may be interpreted as an indirect description of the distribution of slow strain wave velocities in the seismic region.

In a seismic region with non-accelerating cumulative Benioff strain, random matrix model (24) may describe the slow strain wave velocity distribution with a constant value of $\beta$. However, in a seismic activation region, it is necessary to allow for $\tau$ dependence of the damage parameter $\beta \left(\tau \right)$ and function$F(\omega ,\beta (\tau \left)\right)$ to describe the distribution of bound state eigenvalues solving nonlinear eigenvalue problem (19). The renormalization group theory of earthquakes suggests this time dependence can be precisely quantified using statistical mechanics critical scaling theory. Requiring that the critical scaling theory describe statistics of solutions to a nonlinear eigenvalue problem suggests the relevant statistical mechanical system is 2D, and requiring that the bound state eigenvalue statistics account for the Gutenberg-Richter law suggests the relevant statistical mechanical system contains vortex-antivortex dipole lengths R distributed according to an inverse power law at thermal equilibrium.

3. Results

3.1. FDM Random Matrix Eigenvalue Computation

3.1.1. 2D Critical Scaling Theory

In addition to identifying bound state eigenvalues, solutions to Equation (19) also identify scattering resonances of the seismic activation region. These resonant eigenvalues are inherently complex, and with selection of stiffness and damping matrices from a random matrix ensemble, are distributed randomly in the complex plane. For this reason, it is conjectured that a 2D statistical mechanical field theory describes the distribution of resonant eigenvalues in the complex plane.

Furthermore, depending on whether incoming or outgoing boundary conditions are specified on the hemisphere $\mathcal{H}$, resonant eigenvalues may be labeled as positive or negative. This suggests the 2D statistical mechanical system is a 2D Dyson gas, whereby the statistics of bound state eigenvalues describes the size distribution of positive-negative dipoles in the complex plane. From this point of view, a correspondence between metastable cluster and dipole size distribution statistics is established, since a Gutenberg-Richter distribution of earthquake corner frequencies, is:

$$f\left(\omega \right)d\omega \propto {\omega }^{2b-1}d\omega \propto R^{-2b-1}dR,$$

(25)

where $R=1/\omega$, and the distribution of dipole sizes satisfies a power law [26]:

$$f_{dipole}\left(R\right)\propto R^{-2\pi \beta \left(\tau \right)}.$$

(26)

In absence of further development, the renormalization group theory of earthquaes postulates the existence of $\tau$-dependent statistical mechanical systems related by renormalization group transformation, whereby the cumulative Benioff strain at time $\tau$ equates to the expectation value ${\langle C\rangle }_{\tau }$ of a statistical mechanical observable C for which:

$${\displaystyle \frac{d^k{\langle C\rangle }_{\tau }}{d{\tau }^k}},$$

(27)

is infinite for some finite positive integer k at a sequence of values ${\tau }_i$ approaching ${\tau }_0$, which locates a statistical mechanical critical point [13,25]. With stochastic modeling of the seismic activation region presented in this article, the renormalization group theory of earthquakes is developed by conjecturing that the renormalized statistical mechanical system is a 2D Dyson gas, and the renormalization group fixed point is a critical point of the KT phase transition.

A more precise characterization of the 2D statistical mechanical system can be conjectured by noting that if Equation (20) is altered to include $\omega$ dependence of seismic activation region stiffness and damping it becomes:

$$-i\omega \left[\begin{array}{c}\mathbf{u}\\ -i\omega \tilde{M}\mathbf{u}\end{array}\right]=\left[\begin{array}{cc}0& {\tilde{M}}^{-1}\\ -\tilde{K}(\omega ,\tau )& -\tilde{D}(\omega ,\tau ){\tilde{M}}^{-1}\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ -i\omega \tilde{M}\mathbf{u}\end{array}\right],$$

(28)

when $\mathbf{F}(\omega ,\tau )=0$. This equation implies the matrix:

$$\tilde{L}(\omega ,\tau )-\eta I=\left[\begin{array}{cc}-\eta I& {\tilde{M}}^{-1}\\ -\tilde{K}(\omega ,\tau )& -\tilde{D}(\omega ,\tau ){\tilde{M}}^{-1}-\eta I\end{array}\right],$$

(29)

has determinant 0 at $\eta =-i\omega$. Viewed as a $\tau$ dependent relation between $\omega$ and $\eta$, the condition:

$$det(\tilde{L}(\omega ,\tau )-\eta I)=0,$$

(30)

can be ensemble averaged over the $\tau$ dependent distribution of random stiffness and damping matrices to obtain an average relation between $\omega$ and $\eta$, hereby conjectured to be a Riemann surface spectral curve. From this point of view, if the genus g of the average spectral curve remains constant as $\tau \to {\tau }_0$, the renormalization group flow of the 2D statistical mechanical system may be identifiable as the motion of a point in a Siegel moduli space [8]. For sake of specificity, it is further conjectured that this 2D statistical mechanical model of accelerating seismic release has $2g$ bosonic fields, where $2g$ is the number of bound state eigenvalues characterizing the metastable cluster of fault material initiating the mainshock event.

3.2. Geophysical Test

From a geophysical testing point of view, if it is true that the real time evolution of a seismic activation region elastic model preceding a mainshock can be quantified in terms of a 2D statistical mechanics model renormalization group flow, expressible as a nonlinear dynamical system of finite phase space dimension, a geophysical signal processing technique known as singular spectrum analysis should apply to determine this phase space dimension [6]. Therefore, it is conjectured that measurements of relative changes in seismic wave velocity between pairs of seismic stations in a seismic region obtained at regular time intervals during a seismic activation series can be input to a time domain multichannel singular spectrum analysis algorithm to output a finite phase space dimension [17]. More specifically, the number of channels of the algorithm equates to the number of station pairs, and the number of singular values output by the algorithm in different time windows preceding a mainshock event counts the number of unstable stress/strain modes contributing to mainshock rupture nucleation. 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 [24].

4. Discussion

Previous work has identified predicting the time of occurrence of mainshock events as an application of statistical physics models of seismic activation, but this application has not yet been realized [5]. In more recent times, the artificial intelligence algorithm QuakeGPT has been developed for the purpose of forecasting earthquake occurrence, using seismic event record training data created with a stochastic simulator [23]. Therefore, a practical application of statistical mechanics models of seismic activation may be to be improve stochastic simulation of seismic event records for use in earthquake forecasting technology, acknowledging that rigorous tests of model validity against real seismic data must be passed before achieving this objective can be considered a realistic possibility.

In conclusion, work towards improving current earthquake early warning systems can proceed in two directions. Firstly, work can be done to determine whether or not observed changes of the Earth’s elastic velocity model preceding mainshock events can be processed to extract an integer identifiable as the phase space dimension of a nonlinear dynamical system. Secondly, theoretical work can be done to elaborate upon the statistical mechanics model of seismic activation presented in this article to determine other tests of its scientific validity and potential for practical application.