Simulation study of earthquakes based on the two-dimensional Burridge-Knopoff model with the long-range interaction

Spatiotemporal correlations of the two-dimensional spring-block (Burridge-Knopoff) models of earthquakes with the long-range inter-block interactions are extensively studied by means of numerical computer simulations. The long-range interaction derived from an elasticd theory, which takes account of the effect of the elastic body adjacent to the fault plane, falls off with distance r as 1/r^3. Comparison is made with the properties of the corresponding short-range models studied earlier. Seismic spatiotemporal correlations of the long-range models generally tend to be weaker than those of the short-range models. The magnitude distribution exhibits a near-critical'' behavior, i.e., a power-law-like behavior close to the Gutenberg-Richter law, for a wide parameter range with its B-value, B\simeq 0.55, insensitive to the model parameters, in sharp contrast to that of the 2D short-range model and those of the 1D short-range and long-range models where such a near-critical’’ behavior is realized only by fine-tuning the model parameters. In contrast to the short-range case, the mean stress-drop at a seismic event of the long-range model is nearly independent of its magnitude, consistently with the observation. Large events often accompany foreshocks together with a doughnut-like quiescence as their precursors, while they hardly accompany aftershocks with almost negligible seismic correlations observed after the mainshock.

💡 Research Summary

This paper investigates the statistical and spatiotemporal properties of a two‑dimensional Burridge‑Knopoff (BK) spring‑block model when the inter‑block interaction is extended from nearest‑neighbor to a long‑range form derived from elastic theory. By assuming an infinite, homogeneous, isotropic elastic body surrounding the fault plane and using a static approximation (elastic wave speed much larger than rupture speed), the authors obtain an effective spring constant that decays as 1/r³ with the distance r between blocks. The resulting equation of motion (dimensionless) includes a sum over all blocks, each contributing a term proportional to (u_j – u_i)/r⁵, where u denotes block displacement. For comparison, the standard short‑range (nearest‑neighbor) model is recovered when the interaction range is truncated.

The model employs a velocity‑weakening friction law originally introduced by Carlson and Langer: φ(ṽ) = 1 – σ for positive slip velocity, decreasing further with a rate α as velocity increases, while backward slip is prohibited (φ = –∞ for ṽ < 0). The two friction parameters, σ (instantaneous drop) and α (rate of weakening), together with the elastic stiffness ratio l = √(k_c/k_p), control the dynamics. Loading is taken to be infinitesimally slow (ν → 0) so that the system evolves only during slip events, mimicking the quasi‑static tectonic loading of real faults.

Numerical simulations are performed using a fourth‑order Runge‑Kutta integrator with Δt = 10⁻³, on lattices up to 160 × 80 (periodic in the strike direction, free in depth). After discarding an initial transient of 10⁵ events, 10⁵–10⁷ earthquakes are generated for each parameter set, providing robust statistics. The magnitude µ of an event is defined as µ = ln M, where M = Σ_i Δu_i is the total slip summed over all participating blocks.

Key findings include:

1. Magnitude distribution – For a broad range of α (0 ≲ α ≲ ∞) the event‑size distribution R(µ) exhibits a near‑critical power‑law segment with a Gutenberg‑Richter (GR) exponent B ≈ 0.55 (corresponding to b ≈ 0.85). This “near‑critical” behavior persists without fine‑tuning of l, σ, or α, in stark contrast to the short‑range 2D model where B varies strongly with parameters and only a narrow region yields GR‑like scaling. Small α (< 1) suppresses large events, while larger α expands the scaling range.

2. Stress drop – The mean stress drop Δσ̄ is essentially independent of event magnitude, reproducing the observed magnitude‑independent stress release in real earthquakes. In the short‑range model, Δσ̄ grows with magnitude, highlighting the role of long‑range coupling in homogenizing stress redistribution.

3. Spatiotemporal correlations – Large events are frequently preceded by foreshocks and a doughnut‑shaped quiescence zone surrounding the future epicenter. After the mainshock, aftershocks are almost absent, and the seismic time‑correlation function decays rapidly to zero. This pattern mirrors observations of foreshock clustering and limited aftershock activity in some subduction zones, and suggests that long‑range elastic coupling quickly equilibrates the stress field after rupture.

4. Dimensional effects – An appendix treats the one‑dimensional version with a 1/r² interaction (derived by imposing rigidity in the depth direction). The 1D long‑range model still requires parameter fine‑tuning to achieve GR scaling, whereas the 2D long‑range model naturally exhibits robust scaling, underscoring the combined influence of spatial dimensionality and interaction range on universality classes, analogous to equilibrium phase transitions.

Overall, the study demonstrates that incorporating a physically motivated 1/r³ long‑range interaction into the BK framework yields statistical features—GR‑like magnitude distribution, magnitude‑independent stress drop, and realistic foreshock/aftershock patterns—that align closely with empirical seismicity, while reducing the need for delicate parameter adjustment. The work thus advances the realism of spring‑block earthquake models and provides a bridge between simple mechanical analogues and the complex, scale‑invariant behavior observed in nature. Future extensions could include dynamic elastic wave propagation, anisotropic elastic constants, and direct quantitative comparison with regional seismic catalogs to further validate the model’s predictive power.