Preprint
Article

This version is not peer-reviewed.

Seismic Activation Modeling with Statistical Physics

Submitted:

25 July 2025

Posted:

25 July 2025

Read the latest preprint version here

Abstract
A correspondence between fault damage mechanics and critical point models of seismic activation is presented, and a method of testing the correspondence against seismic measurements is outlined.
Keywords: 
;  ;  ;  

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 N 26 N × 119 E 123 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 τ = τ 0 . In this figure, τ is a real time parameter, and τ 0 is the characteristic time of major earthquake recurrence assuming an earthquake of similar magnitude occurred in the same region at τ = 0 [21,29]. 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 τ 0 .
As a means of predicting the time τ = τ 0 at which a mainshock event preceded by seismic activation occurs, it has been hypothesized that the average seismic moment M ¯ τ of earthquakes occuring in time intervals ( τ , τ + Δ τ ) preceding a mainshock event obeys an inverse power of remaining time to failure law:
M ¯ τ 1 ( τ 0 τ ) γ 1
and that the cumulative Benioff strain C ( τ ) , defined as:
C ( τ ) = i = 1 n ( τ ) M 0 , i 1 / 2 ,
where M 0 , i is the seismic moment of the i t h earthquake in the region starting from a time τ = 0 preceding the mainshock event, and n ( τ ) is the number of earthquakes occurring in the region up to time τ , satisfies [27]:
C ( τ ) = a b ( τ 0 τ ) γ 2 , γ 2 = 1 γ 1 / 2 .
The exponent selection of 1 / 2 in equation (2) is not necessary to derive formula (3) with a different arithmetic relation between γ 1 and γ 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 γ 2 is 0.3 [5,28]. 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,31].
A model of seismic activation based on fault damage mechanics (FDM) has been used to derive equation (3) with a value γ 2 = 1 / 3 [2]. 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 τ τ 0 . The result γ 2 = 1 / 3 is obtained from an effective shear modulus evolution equation based on non-equilibrium thermodynamic considerations [1].
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 γ 2 = 1 / 4 [21]. In this derivation, the inverse power of remaining time to failure law:
M ¯ τ 1 ( τ 0 τ ) 3 / 2
is asserted based on identifying the mean rupture length R ¯ τ of earthquakes occuring at time τ with the correlation length of a statistical mechanical system described by a time dependent Ginzburg-Landau equation for which:
R ¯ τ 1 ( τ 0 τ ) 1 / 2 ,
which implies relation (4) given the scaling relation M ¯ τ R ¯ τ 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 with provision of computational evidence. The second objective is to outline how the conjectured correspondence can be tested against seismic measurements.
To meet these objectives, the article develops two different seismicity models: the earthquake cascade model and the renormalization group theory of earthquakes [25,26]. 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 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 2D 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 by finite element analysis. 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 applications of seismic activation statistical physics models.

Materials and Methods

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 ( τ ) as the number of metastable clusters of fault material of area 2 n at time τ , and provides time evolution equations for these numbers based on cluster joining and earthquake occurrence rates. Depending on selection of these rates, either of the earthquake magnitude cumulative distribution curves shown in Figure 1 can be approximated by the model. For example, at sufficiently low rates of earthquake occurrence, time invariant metastable cluster size statistics are:
N 2 = N 1 2 b , N 3 = N 2 2 b , . . . , N n = N n 1 2 b ,
for some constant b , whereby:
N k = 2 k N 1 ( 2 k ) b ,
and:
2 k N k ( 2 k ) 1 b .
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 if b = b 0.5 . To explain this statement, note that if N c M is the number of earthquakes of moment magnitude greater than or equal to M occurring during the time interval ( τ , τ + Δ τ ) , the Gutenberg-Richter law implies:
N c M 10 b M ,
which using the seismic moment-magnitude relation:
M = log 10 M 9 1.5 ,
and corner frequency scaling relation:
M ω 3 ,
implies the frequency f ( ω ) of earthquakes with corner frequency ω and rupture area proportional to ω 2 satisfies:
f ( ω ) ω 2 b 1 = ( ω 2 ) 0.5 b .
To modify the earthquake cascade model so it accounts for a continuum of possible metastable cluster sizes, the approach presented here is to adapt previously developed stochastic process models of landslide occurrence to model earthquake occurrence. Importantly, these stochastic process models depend on whether the landslide is in a condition of secondary or tertiary creep, so their adaptations to describing seismicity depend on whether cumulative Benioff strain is increasing at a constant rate or accelerating.

FDM Random Matrix Model

A scalar Kesten process is a stochastic process:
v k + 1 = λ k + 1 v k + ϵ k + 1 ,
where i is a discrete time index, v 0 = 1 , and λ k and ϵ k are identical and independently distributed random variables. This process has been used to model sequences of landslide velocity measurements by assuming λ k and ϵ k are log-normally distributed with log λ k < 0 , in which case the sequence of velocities are distributed according to an inverse gamma distribution P ( v ) d v satisfying:
P ( v ) v b 1 ,
for some constant b for sufficiently large velocities v [16]. In developing 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 of interest, and output a random variable whose distribution accounts for the Gutenberg-Richter law.
To this end, first suppose the seismic activation region is contained within a hemisphere H centered on the Earth’s surface, and subsurface material within the hemisphere is meshed with finitely many structural elements. Then, according to classical mechanics, a vector of nodal displacements U ( t ) to a finite strain elastoplastic deformation of material within the hemisphere at time τ + t satisfies a differential equation:
M ˜ U ¨ ( t ) + D ˜ ( τ ) U ˙ ( t ) + K ˜ ( τ ) U ( t ) = F ( t ) ,
where M ˜ , D ˜ ( τ ) , and K ˜ ( τ ) are the finite element mass, damping, and tangent stiffness matrices at time τ , F ( t ) is the vector of external forces acting on the nodes, and the time increment t is sufficiently small that the matrix coefficients can be regarded as constant over the time interval ( τ , τ + t ) [3]. Equivalently, differential equation (15) can be written:
d d t U ( t ) M ˜ U ˙ ( t ) = 0 M ˜ 1 K ˜ ( τ ) D ˜ ( τ ) M ˜ 1 U ( t ) M ˜ U ˙ ( t ) + 0 F ( t ) ,
which in absence of external forcing has the solution:
U ( t ) M ˜ U ˙ ( t ) = e L ˜ ( τ ) t U ( 0 ) M ˜ U ˙ ( 0 ) ,
where:
L ˜ ( τ ) = 0 M ˜ 1 K ˜ ( τ ) D ˜ ( τ ) M ˜ 1 .
A special case of solution (17) occurs when e L ˜ ( τ ) t is a symplectic matrix, in which case the equation may be interpreted as a map between tangent spaces of a classical phase space manifold induced by an energy preserving flow. In this case, each eigenvector of L ˜ ( τ ) with eigenvalue 0 identifies a marginally stable nodal displacement of the activation region, in that it is on the cusp of stable oscillation and unstable exponential growth. In accordance with previous fault dynamic investigations, it is these marginally stable states of nodal displacement that constitute a classical mechanical definition of metastable clusters [7]. Assuming activation region fault material is being strained at a non-zero rate, the momentum of the marginally stable states constitutes a classical mechanical definition of metastable cluster size.
Next, to describe the metastable cluster distribution with a stochastic process, define a stochastic process P k of vector random variables to replace the classical mechanistic vectors M ˜ U ˙ ( k Δ t ) , and suppose:
P k + 1 = J k + 1 P k + E k + 1 ,
where J k is an m × m random matrix and E k is a 1 × m random vector, for m equal to the number of marginally stable states. For example, in the special case where J k is diagonal with independent log-normally distributed entries, and E k has independent log-normally distributed entries, process (19) is equivalent to m independent Kesten processes, and these processes determine probability distributions of the momentum magnitudes p k of the marginally stable states whose mean values, in absence of accelerated seismic release, remain constant in time.
Here, instead of supposing the p k are independently distributed, we conjecture stochastic process (19) can be defined to return real positive values p 1 , p 2 , . . . , p m with a random matrix eigenvalue distribution:
P ( p 1 , p 2 , . . . , p m ) = e k = 1 m ln F ( p k , β ( τ ) ) | Δ ( p 1 , p 2 , . . . , p m ) | β ( τ ) ,
where Δ ( p 1 , p 2 , . . . , p m ) is the Vandermonde determinant of the positive eigenvalues, ln F ( p , β ( τ ) ) is a confining potential of the positive eigenvalues in 1 spatial dimension, and β ( τ ) is a time dependent damage parameter [12].

Results

FDM Random Matrix Eigenvalue Computation

2D Critical Scaling Theory

The renormalization group theory of earthquakes postulates the existence of τ -dependent statistical mechanical models related by renormalization group transformation, whereby the cumulative Benioff strain at time τ equates to the expectation value C τ of a statistical mechanical observable C for which:
d k C τ d τ k ,
is infinite for some finite positive integer k = k i at a sequence of values τ i approaching τ 0 , and τ 0 identifies a model fixed by the renormalization group transformation [25]. In the thermodynamic phase diagram of the statistical mechanical model, each value τ i locates a point of phase transition, while the value τ 0 locates a critical point [13].
The renormalization group theory of earthquakes is hereby developed by conjecturing that there exists a sequence of values p ¯ i for which:
d k i ln F ( p ¯ i , β ( τ ) ) d τ k i ,
is infinite at τ = τ i . More specifically, the function F ( p , β ( τ ) ) can be interpreted as a statistical mechanical partition function with fugacity and inverse temperature parameters p and β ( τ ) . Importantly, this conjectural reinterpretation of the time dependent parameter β ( τ ) in equation (20) as an inverse temperature parameter suggests the underlying statistical mechanical model is a 2D Wick rotatation of a random matrix spectral statistics model in 1+1D [18]. For sake of specificity, it is further conjectured that the 2D statistical mechanical model of accelerating seismic release has 2 m bosonic fields, where m is the number of accelerating marginally stable modes in the activation region. An example of such a statistical field theory is provided by the 2D sine-Gordon model [8].

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 τ rather than the frequency domain [24].

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.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Acknowledgments

Thanks to my family for their support throughout completion of this research. Thanks to Dr. Girish Nivarti, Dr. Evans Boney, and Professor Richard Froese for their willingness to entertain communication regarding the content of the article.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Ben-Zion Y. Collective behavior of earthquakes and faults: Continuum-discrete transitions, progressive evolutionary changes, and different dynamic regimes. Reviews of Geophysics2008,46(4).
  2. Ben-Zion Y, Lyakhovsky V. Accelerated seismic release and related aspects of seismicity patterns on earthquake faults. Earthquake processes: Physical modelling, numerical simulation and data analysis Part II2002, 2385–2412.
  3. Bindel DS. Structured and parameter-dependent eigensolvers for simulation-based design of resonant MEMS. Doctoral dissertation. University of California, Berkeley. 2006.
  4. Bletery Q, Nocquet JM. The precursory phase of large earthquakes. Science2023,381(6655) 297–301. [CrossRef]
  5. Bowman D, Ouillon G, Sammis C, Sornette A, Sornette D. An observational test of the critical earthquake concept. Journal of Geophysical Research: Solid Earth1998,103(B10), 24359–24372. [CrossRef]
  6. Broomhead DS, King GP. Extracting qualitative dynamics from experimental data. Physica D: Nonlinear Phenomena1986,20(2-3), 217–236. [CrossRef]
  7. Carlson JM, Langer JS, Shaw BE. Dynamics of earthquake faults. Reviews of Modern Physics1994,66(2), 657. [CrossRef]
  8. Carpentier D. Renormalization of modular invariant Coulomb gas and sine-Gordon theories, and the quantum Hall flow diagram. Journal of Physics A: Mathematical and General1999,32(21), 3865. [CrossRef]
  9. Chen CC. Accelerating seismicity of moderate-size earthquakes before the 1999 Chi-Chi, Taiwan, earthquake: Testing time-prediction of the self-organizing spinodal model of earthquakes. Geophysical Journal International2003,155(1), F1–5. [CrossRef]
  10. Dou L, Cai W, Cao A, Guo W. Comprehensive early warning of rock burst utilizing microseismic multi-parameter indices. International Journal of Mining Science and Technology2018,28(5), 767-74. [CrossRef]
  11. Fortin J, Stanchits S, Dresen G, Gueguen Y. Acoustic emissions monitoring during inelastic deformation of porous sandstone: comparison of three modes of deformation. Pure and Applied Geophysics2009,166(5), 823–41.
  12. Gautié T, Bouchaud JP, Le Doussal P. Matrix Kesten recursion, inverse-Wishart ensemble and fermions in a Morse potential. Journal of Physics A: Mathematical and Theoretical2021,54(25), 255201. [CrossRef]
  13. Goldenfeld N. Continuous Symmetry. In Lectures on phase transitions and the renormalization group; CRC Press, 2018; pp. 335–350.
  14. Hardebeck JL, Felzer KR, Michael AJ. Improved tests reveal that the accelerating moment release hypothesis is statistically insignificant. J. Geophys. Res.2008,113. [CrossRef]
  15. Ito R, Kaneko Y. Physical mechanism for a temporal decrease of the Gutenberg-Richter b-Value prior to a large earthquake. Journal of Geophysical Research: Solid Earth2023,128(12).
  16. Lei Q, Sornette D. A stochastic dynamical model of slope creep and failure. Geophysical Research Letters2023,50(11). [CrossRef]
  17. Merrill RJ, Bostock MG, Peacock SM, Chapman DS. Optimal multichannel stretch factors for estimating changes in seismic velocity: Application to the 2012 Mw 7.8 Haida Gwaii earthquake. Bulletin of the Seismological Society of America2023,113(3), 1077–1090. [CrossRef]
  18. Peskin ME. An introduction to quantum field theory. CRC press. 2018.
  19. Reches ZE. Mechanisms of slip nucleation during earthquakes. Earth and Planetary Science Letters1999,170(4), 475–86. [CrossRef]
  20. Riccobelli D, Ciarletta P, Vitale G, Maurini C, Truskinovsky L. Elastic instability behind brittle fracture. Physical Review Letters2024,132(24). [CrossRef]
  21. Rundle JB, Klein W, Turcotte DL, Malamud BD. Precursory seismic activation and critical-point phenomena. Microscopic and Macroscopic Simulation: Towards Predictive Modelling of the Earthquake Process2001, 2165–2182.
  22. Rundle JB, Turcotte DL, Shcherbakov R, Klein W, Sammis C. Statistical physics approach to understanding the multiscale dynamics of earthquake fault systems. Reviews of Geophysics2003,41(4). [CrossRef]
  23. Rundle JB, Fox G, Donnellan A, Ludwig IG. Nowcasting earthquakes with QuakeGPT: Methods and first results. In Scientific Investigation of Continental Earthquakes and Relevant Studies; Springer Nature, Singapore, 2025.
  24. Sacchi M. FX singular spectrum analysis. In Cspg Cseg Cwls Convention, 392–395.
  25. Saleur H, Sammis C, Sornette D. Renormalization group theory of earthquakes. Nonlinear Processes in Geophysics1996,3(2), 102–109. [CrossRef]
  26. Turcotte DL, Malamud BD. Earthquakes as a complex system. International Geophysics2002,81, Academic Press.
  27. Tzanis A, Vallianatos F, Makropoulos K. Seismic and electrical precursors to the 17-1-1983, M7 Kefallinia earthquake, Greece: Signatures of a SOC system. Physics and Chemistry of the Earth, Part A: Solid Earth and Geodesy2000,25(3), 281–7. [CrossRef]
  28. Vallianatos F, Chatzopoulos G. A complexity view into the physics of the accelerating seismic release hypothesis: Theoretical principles. Entropy2018,20(10), 754.
  29. Vallianatos F, Sammonds P. Evidence of non-extensive statistical physics of the lithospheric instability approaching the 2004 Sumatran–Andaman and 2011 Honshu mega-earthquakes. Tectonophysics2004,590, 52–8. [CrossRef]
  30. Xue L, Belytschko T. Fast methods for determining instabilities of elastic–plastic damage models through closed-form expressions. International Journal for Numerical Methods in Engineering2010,84(12), 1490–518. [CrossRef]
  31. Zhou S, Johnston S, Robinson R, Vere-Jones D. Tests of the precursory accelerating moment release model using a synthetic seismicity model for Wellington, New Zealand. J. Geophys. Res.2006,111. [CrossRef]
Figure 1. Plot of the cumulative distribution of earthquakes of different moment magnitudes in a seismic activation region in two different time intervals of equal width preceding occurrence of a major earthquake at Δ τ = τ 0 τ = 0 [21,29].
Figure 1. Plot of the cumulative distribution of earthquakes of different moment magnitudes in a seismic activation region in two different time intervals of equal width preceding occurrence of a major earthquake at Δ τ = τ 0 τ = 0 [21,29].
Preprints 169670 g001
Table 1. Approximate relation between earthquake magnitude, fault material displacement, and fault rupture length.
Table 1. Approximate relation between earthquake magnitude, fault material displacement, and fault rupture length.
Moment Magnitude Average Fault Material Displacement (m) Fault Rupture Length (km)
4 0.05 1
5 0.15 3
6 0.5 10
7 1.5 30
8 5 100
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.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

Disclaimer

Terms of Use

Privacy Policy

Privacy Settings

© 2026 MDPI (Basel, Switzerland) unless otherwise stated