A Hierarchical Fractional N-body Framework for Coulomb Stress Evolution in Fault Networks: Exact Analytical Solutions and Falsifiable Seismicity Scaling Laws

Chishtie, Farrukh A.

doi:10.5194/egusphere-2026-3498

Preprints

https://doi.org/10.5194/egusphere-2026-3498

Preprints

21 Jul 2026

| 21 Jul 2026

A Hierarchical Fractional N-body Framework for Coulomb Stress Evolution in Fault Networks: Exact Analytical Solutions and Falsifiable Seismicity Scaling Laws

Farrukh A. Chishtie

Abstract. We develop an exact analytical framework for the hierarchical fractional dynamics of interacting many-body systems and apply it to the nonlinear Coulomb stress transfer that governs seismicity in fractally organised fault networks. The framework rests on a scaling relation α_k= 2 − 2 / (N_k+ 1) that connects the order of a Riemann–Liouville fractional evolution equation at hierarchical level k to the number N_k of interacting bodies at that level, derived in a companion paper (Chishtie, 2026, Physics Open) and applied here to the seismogenic setting. Closed-form solutions are obtained via parametric trigonometric representations σ_k(θ) ∝ sin⁴θ at each level together with Chebyshev polynomial inversions of the time–parameter relation, and they converge to the classical wave equation as N_k → ∞. We embed the framework into the Time-Dependent Stress Response (TDSR) seismicity model of Dahm and Hainzl (2022, J. Geophys. Res.) to describe the hierarchical Coulomb stress cascade through a fractally organised fault network. Three quantitative, falsifiable scaling laws follow from the closed-form solutions with no free parameters: an Omori–Utsu aftershock decay exponent p_k = α_k/ 2 = 1 − 1 / (N_k+ 1) that stratifies across generations; a spatial Coulomb stress falloff ℓ^−(1+α_k)that departs measurably from the classical elastic ℓ⁻³law; and a Gutenberg–Richter b-value b_k ≈D_f N_k/ (N_k+ 1) that accounts for the longstanding discrepancy between the fractal-dimension expectation b = D_f≈ 1.5 – 2 and the commonly observed b ≈ 1 as a finite-N_k fractional correction. Numerical verification using the Grünwald–Letnikov scheme against the Mittag-Leffler exact series solution of the fractional relaxation equation confirms the analytical results. The standard TDSR model is recovered exactly as N_k → ∞ at all hierarchical levels and is therefore a special case of the present framework. We present illustrative qualitative comparisons with the 2023 Kahramanmaraş doublet and the 2024 Noto Peninsula earthquake; a systematic generation-stratified analysis across many sequences is set out as the primary observational test of the framework.

Received: 15 Jun 2026 – Discussion started: 21 Jul 2026

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

Download & links

Farrukh A. Chishtie

Status: final response (author comments only)

RC1:
'Comment on egusphere-2026-3498', Anonymous Referee #1, 17 Aug 2026

The manuscript requires a major revision for the following reasons:
1-The introduction requires a revision. Motivations and main objectives of the manuscript are not clear enough.
2-Definition 1 requires a rewritten.
3-Equation 10 requires a revision. Check dimensions please.
4-Equation 10 is not new. It has been treated largely in the literature.
5-Equation 11 is weird enough. Please check very carefully.
6-The confrontations with observations with observations.
7-Discuss the implications of fractional calculus in earthquakes. Check many missing references:
https://doi.org/10.3390/geohazards6030052; https://www.aimspress.com/article/doi/10.3934/math.2025119; https://doi.org/10.1016/j.soildyn.2018.09.006; https://doi.org/10.3390/fractalfract9050320; https://doi.org/10.1177/1077546314557553; https://doi.org/10.1007/s00707-021-03128-9 https://doi.org/10.1007/s00707-022-03213-7; https://doi.org/10.1016/j.strusafe.2024.102525; https://doi.org/10.1088/1742-2140/aabe61; https://doi.org/10.1080/01495739.2021.1919585; https://doi.org/10.1007/s11803-002-0070-5; https://doi.org/10.1007/s00707-020-02929-8; https://doi.org/10.1016/j.soildyn.2018.09.006;
Please revise.

Citation: https://doi.org/10.5194/egusphere-2026-3498-RC1
- AC1: 'Reply on RC1', Farrukh Chishtie, 14 Sep 2026
  
  I thank the referee for a careful and constructive report. I accept the recommendation of major revision. In preparing this response I re-examined the governing equation of Section 3 and its link to the seismicity rate numerically, and comments 3–5 led me to two substantive corrections, described under those points. The changes I will make are as follows.
  1. Introduction. I will rewrite the Introduction so that the motivation and objectives are stated in one explicit paragraph. The motivation is that existing physics-based seismicity models (Coulomb failure, rate-and-state, and their TDSR synthesis) treat the fault population as an effective medium and give no closed-form description of stress transfer through a hierarchically organised fault network, while existing fractional extensions fit the fractional order to data. The objectives are (i) to fix the orders of the fractional stress-transfer operator from the number of interacting fault segments at each hierarchical level, using the selection rule derived in the companion paper; (ii) to embed the resulting stress evolution in the TDSR rate framework of Dahm and Hainzl (2022); and (iii) to derive from it spatial and temporal scaling laws that can be falsified by generation-stratified catalog analysis.
  2. Definition 1. I will rewrite the definition of the Riemann–Liouville operator with the function space, the range n − 1 < α ≤ n, the associated fractional initial-value structure, and its relation to the Caputo operator stated in full, and I will add a table of all symbols used in Sections 2–4. I will state explicitly which initial-value data are physically prescribed for a fault system with pre-existing deformation.
  3. Equation (10): dimensions. The referee is correct. As written, the generalised diffusivity must carry units of m² s^(−α_k) and the source term units of Pa s^(−α_k), and neither is stated. In the revision every coefficient will be given with its units, and the equation will be written with an explicit level-k length scale L_k and time scale τ_k so that all fractional exponents act on dimensionless ratios.
  4. Equation (10): novelty. I agree. The time-fractional diffusion-wave equation is due to Schneider and Wyss (1989) and Mainardi (1996) and is reviewed by Metzler and Klafter (2000) and Mainardi, Luchko and Pagnini (2001). The revised text will cite this literature and describe the equation as adopted rather than derived. The contribution of the paper is that the fractional orders are fixed by the integer N_k rather than fitted, and the text will be revised so that this is the stated claim.
  5. Equation (11), and two corrections to the formulation. The referee's instinct is right, and the problem is deeper than the mislabelled inversion. Equation (11) is dimensionally inconsistent (a_k²/D_k^(1/2) is not a time) and gives the forward map t(θ), not the inversion. In checking it I found two further problems.
  First, the spatial scaling law of Eq. (16), Δσ ∝ ℓ^(−(1+α_k)), does not follow from Eq. (10) as written: a time-fractional equation with an integer Laplacian has a stretched-exponential far-field profile, not a power law. A power-law tail ℓ^(−(d+α)) is the signature of a Riesz operator of order α in space. This is the operator of the companion paper, in which α_k = 2 − 2/(N_k+1) is the spatial (Riesz) order and the conjugate temporal order is μ_k = 2/(N_k+1), with α_k + μ_k = 2. In the submitted manuscript α_k was placed on the time derivative, which is the origin of the inconsistency the referee detected. Section 3 will be rewritten around the space–time fractional equation
  ₀^RL D_t^(μ_k) σ_c^(k) = −D_k (−Δ)^(α_k/2) σ_c^(k) + S_k, μ_k = 2/(N_k+1), α_k = 2 − 2/(N_k+1),
  with the Riemann–Liouville initial condition ₀D_t^(μ_k−1) σ_c^(k)|₀₊ = φ^(k) prescribed by the pre-existing stress state, and with the solution and its asymptotics obtained from the known fundamental solution (Mainardi, Luchko and Pagnini, 2001) rather than asserted. Equation (11) will be removed; the parametric trigonometric solution of the companion paper will be presented as the level-k inter-segment separation only, and its role in setting the initial data will be made explicit. I have verified numerically, by Mittag-Leffler evaluation of the Fourier-inverted fundamental solution and an independent FFT check of the spatial profile, that under this formulation the far-field stress decays as ℓ^(−(1+α_k)) with the predicted amplitude. Eq. (16) is therefore retained and now follows from the governing equation.
  Second, the submitted manuscript identified the Omori exponent with the decay exponent of the Coulomb stress perturbation. In re-examining that step against the TDSR rate integral (14), I find it is not supported: for both a uniform and a steady-state source density, a stress step followed by power-law decay yields an asymptotic Omori exponent of unity independent of the decay exponent, with the decay affecting only the transient; for the steady-state case this reproduces the classical result of Dieterich (1994). The revised manuscript will therefore withdraw p_k = 1 − 1/(N_k+1) as a derived result. In its place, the corrected governing equation yields the aftershock migration exponent, L(t) ∝ t^(μ_k/α_k) = t^(1/N_k), which is measurable directly from relocated catalogs and will be presented as the primary temporal scaling law alongside the spatial law. The Omori exponent will be treated as an output of the TDSR rate that depends on the stress-distance distribution of the level-k sources; the framework's implication for that distribution will be stated as a hypothesis at the same level as the b-value relation and tested as such. Table 2, Figs. 2–5, Table 1 and the abstract will be revised accordingly.
  6. Confrontation with observations. I read this comment as asking for a quantitative rather than illustrative comparison. Section 5 will be replaced by generation-stratified estimates of the migration exponent, the spatial falloff, and the b-value, with uncertainties, for at least one relocated sequence using nearest-neighbour cluster families, and will state clearly which predictions are supported, which are not, and which remain untested. The predictions are fixed before that analysis is performed and will not be adjusted to fit it.
  7. Fractional calculus in earthquake science and missing references. I will add a subsection surveying fractional-calculus applications in earthquake science that distinguishes the fractional viscoelastic and structural-response literature the referee lists (fractional response spectra, fractional single-degree-of-freedom damping, base isolation, thermoelastic and structural-safety formulations) from the fractional treatment of seismicity and stress transfer to which the present paper belongs, and will cite the works listed. I note that one DOI in the list is repeated and that the GeoHazards paper is a statistical analysis of the 2025 Tokara swarm rather than a fractional-calculus study; I will cite it where swarm statistics are discussed.
  I am grateful to the referee. Comments 3–5 identified a genuine inconsistency, and the corrected formulation is both more constrained and more consistent with the companion paper than the one submitted.
  
  Citation: https://doi.org/10.5194/egusphere-2026-3498-AC1
RC2:
'Comment on egusphere-2026-3498', Anonymous Referee #2, 09 Sep 2026

My comments on the manuscript, entitled ‘A Hierarchical Fractional N-body Framework for Coulomb Stress Evolution in Fault Networks: Exact Analytical Solutions and Falsifiable Seismicity Scaling Laws’ by Farrukh A. Chishtie

(MS No.: egusphere-2026-3498)

Although the study is interesting and significant, I cannot follow the manuscript because the so-called N-body model was not described and numerous important terms were not clearly defined and explained. I suggested that the author should revise the manuscript and then re-submit it.

Problems:
(1)The author should plot a diagram to display his N-body model.

(2)On Line 48: The author should clearly explain the definitions of S_o and the skin parameter ds.

(3)On Line 50: What is the definition of χ(ζ, t), i.e., the time-evolving density of fault sources at stress level ζ.

(4)On Line 77: What is the parameter q?

(5)What is the relation between the N-body model and the seismogenic Volume mentioned on Line 99?

(6)Could the author clearly describe the meaning of ‘level’ in sub-section: 2.1 Geometry and definitions?

Citation: https://doi.org/10.5194/egusphere-2026-3498-RC2
- AC2: 'Reply on RC2', Farrukh Chishtie, 14 Sep 2026
  
  I thank the referee for the report and agree that the model must be described on the page rather than by reference to the companion paper. I note at the outset that, in responding to Referee #1, I have found it necessary to reformulate the governing equation of Section 3 (see AC1, point 5); the revisions below are written for the reformulated version, and several of the referee's questions are easier to answer clearly under it.
  (1) Diagram of the N-body model. A new Figure 1 will show the hierarchical fault network as an N-body system: the seismogenic volume V; the master segment at level 0; the N₁ daughter segments at level 1 and the N₂ sub-segments at level 2 that each of them loads; the effective mass, area and separation assigned to each segment; and the direction of Coulomb stress transfer from one level to the next. A second panel will show how the stress perturbation from a level-(k−1) rupture enters the level-k subsystem as the source term of the fractional stress equation.
  (2) S₀ and δσ (line 48). S₀ is the Coulomb stress acting on a fault element and δσ is the stress-skin parameter of the TDSR model (Dahm and Hainzl, 2022): the stress scale over which the mean time to failure changes by a factor e, so that ζ = S₀ − σ_c is the distance to failure in stress units. I suspect these symbols rendered as "So" and "ds" in the referee's viewer; in any case both will be defined in words at first use and listed in a new notation table.
  (3) χ(ζ, t) (line 50). χ(ζ, t) is the number density of fault elements per unit stress distance ζ from failure, so that χ(ζ, t) dζ is the number of elements within dζ of the threshold at time t and its integral over ζ is the total number of elements in V. Elements are advected in ζ by the Coulomb stress history and removed by failure at rate 1/t̄_f(ζ). This will be stated explicitly, together with the initial distribution assumed at each level.
  (4) The parameter "q" (line 77). The symbol at line 77 is θ, the dimensionless parameter of the trigonometric representation σ(θ) ∝ sin⁴θ from the companion paper, which increases monotonically with physical time; it is not a free parameter. In the reformulated Section 3 this representation describes the level-k inter-segment separation and fixes the initial data of the stress equation; it is no longer used as the stress solution itself, and the time–parameter relation of Eq. (11) is removed. θ will be defined at first use.
  (5) Relation between the N-body model and the seismogenic volume V (line 99). V is the domain of the model. The fault elements in V are partitioned by size into levels, and the N_k elements at level k that receive stress from a single level-(k−1) rupture form one N_k-body subsystem. The N-body model is therefore the interaction structure of the elements within V, not a separate object; the number of subsystems at level k, and hence the weight w_k in the rate integral, follows from the branching of the hierarchy. This will be stated in Section 2.1 and shown in the new figure.
  (6) Meaning of "level" (Section 2.1). "Level" is the generation index k in the hierarchy: level 0 is the mainshock segment, level 1 comprises the segments it loads directly, level 2 those loaded by level-1 ruptures, and so on, with characteristic length L_k = L₀ λ^k. I will define this at the start of Section 2.1, use "level" and "generation" consistently, and state that the index labels the geometric position of a fault patch in the hierarchy, not the magnitude of any event.
  I am grateful to the referee; the requested definitions and figure will make the revised manuscript self-contained.
  
  Citation: https://doi.org/10.5194/egusphere-2026-3498-AC2
RC3:
'Comment on egusphere-2026-3498', Anonymous Referee #1, 16 Sep 2026

Missing many references in the field mainly related to scaling law, fractal dimension, and fractional calculus
https://doi.org/10.1016/j.optlastec.2026.116233
https://doi.org/10.3390/app15052759
https://doi.org/10.1016/j.epsl.2026.120225
https://doi.org/10.1016/j.comnet.2026.112649
https://doi.org/10.1007/s00707-020-02929-8
https://doi.org/10.1002/2017JB014927
https://doi.org/10.1007/s00707-021-03128-9

Citation: https://doi.org/10.5194/egusphere-2026-3498-RC3
- AC3: 'Reply on RC3', Farrukh Chishtie, 16 Sep 2026
  
  Thanks for sharing this, I'll include them, plus highlight why my approach moves the fractional calculus forward as well.
  
  Citation: https://doi.org/10.5194/egusphere-2026-3498-AC3
- AC4: 'Reply on RC3', Farrukh Chishtie, 16 Sep 2026
  
  Thanks for sharing this, I'll include them, plus highlight why my approach moves the fractional calculus forward as well.
  
  Citation: https://doi.org/10.5194/egusphere-2026-3498-AC4
- AC5: 'Reply on RC3', Farrukh Chishtie, 16 Sep 2026
  
  Thanks for sharing this, I'll include them, plus highlight why my approach moves the fractional calculus forward as well.
  
  Citation: https://doi.org/10.5194/egusphere-2026-3498-AC5

Farrukh A. Chishtie

Viewed

Total article views: 218 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
142	52	24	218	19	18

HTML: 142
PDF: 52
XML: 24
Total: 218
BibTeX: 19
EndNote: 18

Views and downloads (calculated since 21 Jul 2026)

Month	HTML	PDF	XML	Total
Jul 2026	80	16	10	106
Aug 2026	23	13	2	38
Sep 2026	39	23	12	74

Cumulative views and downloads (calculated since 21 Jul 2026)

Month	HTML	PDF	XML	Total
Jul 2026	80	16	10	106
Aug 2026	23	13	2	38
Sep 2026	39	23	12	74

Viewed (geographical distribution)

Total article views: 206 (including HTML, PDF, and XML) Thereof 206 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 22 Sep 2026

Short summary

Earthquake faults branch into smaller and smaller pieces, like cracks within cracks. We built a mathematical model linking how many fault pieces interact to how quickly aftershocks fade, how far the stress from a quake spreads, and how earthquake sizes are distributed. The model needs no tuned settings and agrees with two recent large earthquakes, in Turkey in 2023 and Japan in 2024. It offers a simpler, physics-based way to understand aftershock patterns and earthquake hazard.


Total:	0
HTML:	0
PDF:	0
XML:	0