Prerpints.org logo
Preprint
Article

This version is not peer-reviewed.

Seismic Activation Modeling with Statistical Physics

Submitted:

26 May 2025

Posted:

27 May 2025

Read the latest preprint version here

Abstract
Starting from a quasi-elastostatic description of a seismic activation region, it is explained how accelerating seismic release preceding a mainshock event can be described in terms of a 2D linear sigma model with quartic self interaction. This model is demonstrated to account for the discrepancy in accelerating release critical exponents calculated by fracture damage mechanics and critical point models of seismic activation.
Keywords: 
seismic activation
;  
fault dynamics
;  
statistical physics
;  
signal processing
Subject: 
Physical Sciences  -   Applied Physics

Introduction

An increase in the number of intermediate sized earthquakes ( M > 3 . 5 ) in a seismic region preceding the occurrence of an earthquake with magnitude M > 6 , referred to as seismic activation, has been documented by various researchers [6]. 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 [11]. 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 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 τ 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 intervals of time ( τ , τ + Δ τ ) 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 has 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 [6,28]. Notably, validity of equation (1) has been questioned by some researchers who claim measurements of seismic activation can be explained in terms of main event foreshock and aftershock occurrence without acceleration of seismic release [16,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 [3]. 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 γ 2 = 1 / 3 is obtained from an equation for time evolution of the shear modulus derived from non-equilibrium thermodynamic considerations [2].
In addition to the FDM model of seismic activation, an empirical statistical physics 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 L ( τ ) of earthquakes occuring at time τ with the correlation length of a statistical physical system described by Ginzburg-Landau mean field theory with a τ -dependent temperature parameter, whereby:
L ( τ ) ∝ 1 ( τ 0 − τ ) 1 / 2 ,
and relation (4) follows from the scaling relation 〈 M 〉 τ ∝ L ( τ ) 3 [22]. Table 1 shows typical fault material displacements and rupture lengths for earthquakes of different moment magnitudes.
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
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 physics formalism [25]. Therefore, the first objective of this article is to clarify how the FDM and CP models of seismic activation can be in correspondence with each other. The second objective of the article is to use this correspondence to advance rigorous testing of seismic activation model predictions against seismic measurements, and in the event of positive experimental verification, advance earthquake prediction technology.
Motivating the presented correspondence between FDM and CP seismic activation models is previous work demonstrating statistical physics renormalization group flow equations can, in certain cases, be identified with differential equations such as the Kolmogorov-Petrovsky-Piskunov (KPP) equation, an equation which has been used to model the time and space dependent distribution of aftershocks in a seismic region following a mainshock event [10,15]. This theoretical work may have application to earthquake prediction if it is true that dimensional reduction of statistical physics models at critical points can be used to systematize dimensional reduction of fault dynamic models in windows of time preceding a mainshock event.
The outline of the article is as follows. Section 2 explains how accelerating seismic release can be described in terms of a 2D statistical physics, and why this description is physically reasonable. Section 3 concludes by commenting on how validity of statistical physics modeling of seismic activation can be tested against seismic measurements.

Materials and Methods

Seismic Activation Fault Dynamics

Figure 2 shows a 2D schematic of earthquake occurrence in a seismic activation region [18]. In this figure, the activation region is shown at 4 different times up to and including the moment after a mainshock event has occurred. At each time, black lines indicate fault ruptures associated with earthquakes that have occurred, and red lines indicate faults where shear strain is accumulating prior to earthquake occurrence. Qualitatively, the picture suggests the occurrence of successively larger earthquakes, associated with successively longer rupture lengths, leads to increased strain along the mainshock fault as seismic activation proceeds. From an FDM point of view, this increased strain occurs with a reduction in the average shear modulus of material in the vicinity of the fault, until fault rupture occurs at time τ = τ 0 .
Quantitatively, this picture of seismic activation leading to rupture along a mainshock fault is supported by modeling of earthquake fault dynamics in 1+1 spacetime dimensions, whereby the differential equation:
A ∂ τ 2 U ( τ , z ) − B ∂ z 2 U ( τ , z ) + C ∂ τ U ( τ , z ) = − sin ( U ( τ , z ) / D ) .
has been used to model both creep along an earthquake fault and rupture propagation, depending on whether or not frictional forces dominate the fault dynamics and shear stress evolution along the fault is more appropriately described with a reaction diffusion equation or a solitary wave equation [8]. In this equation, τ is real time, z coordinates a direction of creep or slip along an earthquake fault, U ( τ , z ) is the local displacement of elastic material across the earthquake fault, A ∂ τ 2 U ( τ , z ) is the local inertial force acting on the fault material, B ∂ z 2 U ( τ , z ) is the local elastic restoring force acting on the fault material, and C ∂ τ U ( τ , z ) and sin ( U ( τ , z ) / D ) are local frictional forces acting on the fault material attributed to contact of the material with tectonic plates on either side of the fault. For C = 0 , an (anti-kink) soliton solution to equation can be interpreted as propagation of earthquake fault rupture [30].
To generalize this description of fault creep and rupture in 1 spatial dimension to 3 spatial dimensions, assume that for τ < τ 0 , material constituting the seismic activation region undergoes a quasi-elastostatic finite strain deformation, whereby at any moment in time it exists in an elastostatic equilibrium configuration in which strain energy is minimized . With this assumption, if the seismic activation region is ascribed a finite element mesh, a nodal displacement ψ → of the region’s equilibrium configuration at tme τ increases the strain energy of the region by:
Δ E = 1 2 ψ → T K s a r ( τ ) ψ → ,
for K s a r ( τ ) equal to the positive definite stiffness matrix of the region at time τ . For τ < τ 0 , K s a r ( τ ) has N eigenvalues λ 0 ( τ ) ≤ λ 1 ( τ ) ≤ … ≤ λ N − 1 ( τ ) , where N is the number of finite element mesh nodes. At τ = τ 0 , K s a r ( τ ) has at least one zero eigenvalue λ 0 identifying a marginally stable seismic displacement ψ → = u → 0 that describes the mainshock faulting mechanism [9]. If we now further suppose that in the limit τ → τ 0 , a subset of the eigenvalues of K s a r ( τ ) , including λ 0 ( τ ) , have eigenvectors with 0 nodal displacement outside of a seismic activation subregion containing the mainshock fault, it follows that the subregion has a stiffness matrix K f a u l t ( τ ) that may be isospectral for τ ≥ τ 0 with real time evoution determined by a solitonic Lax pair [17].

Statistical Physics Critical Scaling Theory

From a classical deterministic view of the seismic activation region, its stiffness matrix K s a r ( τ ) eigenvalues λ 0 ( τ ) , λ 1 ( τ ) , … λ N − 1 ( τ ) undergo a deterministic motion of N points on the real line. However, if we instead view the evolution of the seismic activation region elastic model as a stochastic process in which K s a r ( τ ) is selected from a τ -dependent ensemble of random matrices, a corresponding description of the stiffness matrix eigenvalue evolution must be probabilistic. Therefore, the objective herin is to investigate whether or not real time evolution of the stiffness matrix eigenvalue distribution can be described in terms of statistical physics models, whereby convergence of the stiffness matrix eigenvalue distribution to an isopsectral distribution at τ = τ 0 corresponds to invariance of statistical physics model coefficients at a renormalization group fixed point.
Qualitatively, the eigenvalues λ i ( τ ) of the stiffness matrix K s a r ( τ ) depend on the material composition of the activation region, the level of fracture induced damage of material, and the level of stress applied to material, with fracture induced damage dominating relative elastic parameter changes at seismic wavelengths short in comparison to the mainshock rupture length throughout activation, and applied stress dominating relative elastic parameter changes at long seismic wavelengths just preceding the mainshock. In principle, the real time evolution of the activation region elastic model may be governed by a Boltzman kinetic equation describing how the density of fractures of different lengths in the seismic activation region varies in space and time as fractures grow, propagate, and fuse together.
Quantitatively, if the fracture length distribution in each material element at each point in time is characterized in terms of a set of frequency dependent elastic parameters, such as the shear modulus μ if the activation region material is elastically isotropic, the spatial average of each elastic parameter across the activation region at each time τ specifies a function of frequency [4]. For our purpose, defining P = 1 / ω , and assuming material isotropy and that μ ( P , τ ) is a real number, it is conjectured that μ ( P , τ ) d P = | ϕ ( P , τ ) | 2 d P , where ϕ ( P , τ ) satisfies a Klein-Gordon equation with ϕ 4 interaction in 1+1 ( P , τ ) dimensions describing τ -evolution of a bosonic field false ground state [12]. This conjecture allows for previously reported claims that progression of seismic activation can be described in terms of statistical physics critical scaling theory, assuming the spatial domain of the statistical physics model is 2D, and the statistical physics model is a 2D linear sigma model with quartic self interaction related to the 1+1D bosonic field theory by Wick rotation [20].

Results

The 2D linear sigma model with quartic self interaction is described by the Landau free energy functional:
L F = ∫ d 2 r 1 2 ( ∇ ϕ ) 2 + a ( | ϕ | 2 − b ) 2 ,
where the constant b ∝ ( T 0 − T ) is zero at temperature T = T 0 [14]. The mean field Landau free energy density associated with this functional is:
a ( | 〈 ϕ 〉 | 2 − b ) 2 ,
where 〈 ϕ 〉 is the mean field value of ϕ ( r ) . This theory implies mean field scaling:
| 〈 ϕ 〉 | 2 ∝ ( T 0 − T ) ,
that is valid for temperatures far enough from the critical temperature T 0 that the Ginzburg criterion is satisfied.
In terms of elasticity, mean field scaling relation (10) implies:
〈 μ ( P , τ ) 〉 ∝ ( τ 0 − τ ) .
Therefore, assuming 〈 μ ( P , τ ) 〉 is proportional to the average corner frequency of a seismic activation earthquake occurring at time τ , and earthquake corner frequency is inversely proportional to rupture length, the CP rupture length scaling relation:
L ( τ ) ∝ 1 ( τ 0 − τ ) 1 / 2 ,
is obtained. Appreciating that the CP model is a mean field theoretic description of accelerating seismic release, it now follows that the difference between CP and FDM accelerating seismic release critical scaling exponents may derive from anomalous scaling of a 2D statistical physics model underlying the effective 2D linear sigma model that is valid in the critical region T ≈ T 0 where the Ginzburg criterion is violated.

Discussion

Previous research 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 [6]. 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 [5,13,23]. Therefore, a practical application of statistical physics 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.
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 using a 2D statistical physics 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 [7]. More specifically, it is suggested that measurements of relative changes in seismic wave velocity be performed between pairs of seismic stations in a seismic region at regular time intervals during a seismic activation series, and used as input to a time domain multichannel singular spectrum analysis algorithm [19]. The number of channels of this algorithm should equate to the number of station pairs, and the number of singular values output by the algorithm in different time windows preceding a mainshock event should count the number of unstable stress/strain modes contributing to rupture nucleation 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 τ rather than the frequency domain [24].
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, 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.

Author Contributions

Not applicable.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

We encourage all authors of articles published in MDPI journals to share their research data. In this section, please provide details regarding where data supporting reported results can be found, including links to publicly archived datasets analyzed or generated during the study. Where no new data were created, or where data is unavailable due to privacy or ethical restrictions, a statement is still required. Suggested Data Availability Statements are available in section “MDPI Research Data Policies” at https://www.mdpi.com/ethics.

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.

Abbreviations

The following abbreviations are used in this manuscript:
MDPI Multidisciplinary Digital Publishing Institute
DOAJ Directory of open access journals
TLA Three letter acronym
LD Linear dichroism

References

  1. Bena I, Droz M, Lipowski A. Statistical mechanics of equilibrium and nonequilibrium phase transitions: the Yang–Lee formalism. International Journal of Modern Physics B 2005, 19(29), 4269–329.
  2. Ben-Zion Y. Collective behavior of earthquakes and faults: Continuum-discrete transitions, progressive evolutionary changes, and different dynamic regimes. Reviews of Geophysics 2008, 46(4).
  3. 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 II 2002, 2385–2412.
  4. Bindel D.S. Structured and parameter-dependent eigensolvers for simulation-based design of resonant MEMS. Doctoral dissertation. University of California, Berkeley. 2006.
  5. Böse M, Andrews J, Hartog R, Felizardo C. Performance and next-generation development of the finite-fault rupture detector (FinDer) within the United States West Coast ShakeAlert warning system. Bulletin of the Seismological Society of America 2023, 113(2), 648–63.
  6. Bowman D, Ouillon G, Sammis C, Sornette A, Sornette D. An observational test of the critical earthquake concept. Journal of Geophysical Research: Solid Earth 1998, 103(B10), 24359–24372.
  7. Broomhead DS, King GP. Extracting qualitative dynamics from experimental data. Physica D: Nonlinear Phenomena 1986, 20(2-3), 217–236.
  8. Bykov VG. Solitary waves on a crustal fault. Volcanology and Seismology 2001, 22(6), 651–661.
  9. Carlson JM, Langer JS. Properties of earthquakes generated by fault dynamics. Physical Review Letters 1989, 62(22), 2632.
  10. Carpentier D, Le Doussal P. Disordered XY models and Coulomb gases: renormalization via traveling waves. Physical review letters 1998, 81(12), 2558.
  11. 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 International 2003, 155(1), F1–5.
  12. Coleman S. Fate of the false vacuum: Semiclassical theory. Physical Review D 1977, 15(10), 2929.
  13. Cua G., Heaton T. The Virtual Seismologist (VS) method: A Bayesian approach to earthquake early warning. In Earthquake Early Warning Systems, Eds. Gasparini P, Manfredi G, and Zschau J; Springer: Berlin and Heidelberg, Germany, 2007; pp. 97–-132.
  14. Goldenfeld N. Continuous Symmetry. In Lectures on phase transitions and the renormalization group; CRC Press, 2018; pp. 335–350.
  15. Guglielmi AV, Klain BI, Zavyalov AD, Zotov OD. A Phenomenological theory of aftershocks following a large earthquake. Journal of Volcanology and Seismology 2021, 15(6), 373–8.
  16. Hardebeck JL, Felzer KR, Michael AJ. Improved tests reveal that the accelerating moment release hypothesis is statistically insignificant. J. Geophys. Res. 2008, 113.
  17. Kasman A. Glimpses of soliton theory: the algebra and geometry of nonlinear PDEs; American Mathematical Society, 2023.
  18. Lei Q, Sornette D. Anderson localization and reentrant delocalization of tensorial elastic waves in two-dimensional fractured media. Europhys. Letters 2022, 136(3), 1–7.
  19. 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 America 2023, 113(3), 1077–1090.
  20. Peskin ME. An introduction to quantum field theory. CRC press. 2018.
  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 Process 2001, 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 Geophysics 2003, 41(4).
  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 Geophysics 1996, 3(2), 102–109.
  26. Schoenberg M. Elastic wave behavior across linear slip interfaces. The Journal of the Acoustical Society of America 1980, 68(5), 1516–21.
  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 Geodesy 2000, 25(3), 281–7.
  28. Vallianatos F, Chatzopoulos G. A complexity view into the physics of the accelerating seismic release hypothesis: Theoretical principles. Entropy 2018, 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. Tectonophysics 2004, 590, 52–8.
  30. Wu ZL, Chen YT. Solitary wave in a Burridge-Knopoff model with slip-dependent friction as a clue to understanding the mechanism of the self-healing slip pulse in an earthquake rupture process. Nonlinear Processes in Geophysics 1998, 5(3), 121-5.
  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.
Preprints 161038 g001
Figure 1. Plot of the cumulative distribution of earthquakes of different moment magnitudes in a seismic zone 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 zone in two different time intervals of equal width preceding occurrence of a major earthquake at Δ τ = τ 0 − τ = 0 [21,29].
Preprints 161038 g001
Preprints 161038 g002
Figure 2. Schematic illustration of seismic activation in a 2D geometry at four different times τ in which each black line represents an earthquake fault rupture that has already occured, and the red lines represent earthquake faults along which shear stress is increasing prior to rupture [18].
Figure 2. Schematic illustration of seismic activation in a 2D geometry at four different times τ in which each black line represents an earthquake fault rupture that has already occured, and the red lines represent earthquake faults along which shear stress is increasing prior to rupture [18].
Preprints 161038 g002
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.

facebook logotwitter logolinkedin logo
bsky
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

Privacy Settings

© 2025 MDPI (Basel, Switzerland) unless otherwise stated