Beyond static forecasts: a dynamic stress gradient framework for high-resolution aftershock prediction and mitigation

Liao, Boi-Yee

doi:10.5194/npg-33-123-2026

Articles | Volume 33, issue 1

https://doi.org/10.5194/npg-33-123-2026

© Author(s) 2026. This work is distributed under
the Creative Commons Attribution 4.0 License.

https://doi.org/10.5194/npg-33-123-2026

© Author(s) 2026. This work is distributed under
the Creative Commons Attribution 4.0 License.

Articles | Volume 33, issue 1

Research article

|

19 Mar 2026

Research article |

| 19 Mar 2026

Beyond static forecasts: a dynamic stress gradient framework for high-resolution aftershock prediction and mitigation

Boi-Yee Liao

Download

Final revised paper (published on 19 Mar 2026)
Preprint (discussion started on 27 Oct 2025)

Interactive discussion

Status: closed

CC1:
'Comment on egusphere-2025-3348', Jeen-Hwa Wang, 27 Nov 2025

My comments on the manuscript, entitled ‘Beyond Static Forecasts: A Dynamic Stress Gradient Framework for High-Resolution Aftershock Prediction and Mitigation’ (egusphere-2025-3348 (NPG)) authored by Boi-Yee Liao
The author constructed a framework by integrating the rate-and-state friction law, a KPP-type reaction-diffusion equation, and the Banach fixed-point theorem. He applied this framework to study the evolution of Coulomb stress following a large earthquake from numerical simulations of the dynamic redistribution of stress in space and time.
Although the study is interesting and significant, I cannot correctly evaluate whether manuscript may be accepted for publication in its current form or not due to the following basic reasons in ‘Problems for Sciences.’ These reasons make me unable to judge whether the numerical results made by the author are significant for the generation of aftershocks after the 2018 Hualien earthquake in Taiwan or not. Hence, I suggest that the author should answer the questions. This manuscript should be substantially revised and then reviewed again.

Major Problems with Sciences:
(1)In the author’s framework, the rate-and-state friction law, the KPP-type reaction-diffusion equation, and the Banach fixed-point theorem are all one-dimensional. However, his simulations are made on a two-dimensional rectangle. The author should describe the two-dimensional model because the values of model parameters are often distinct in the different axes.
(2)The author wrote: ‘To ensure numerical convergence and model reliability, we adopt spatial grid resolutions of Δx=1320 m and Δy=2600 m, yielding an effective length scale of Δ =2900.’ In numerical simulations, the grid sizes can remarkably influence numerical convergence and stability. The selection of grid sizes for numerical simulations is important. How did the author select the two values: from a theoretical analysis or from numerical tests?
(3)For the study area represented by a rectangle in Figures 5-7, the differences in degrees are 0.5^o (or 55 km) along the longitude and 1^o (or 110 km) along the latitude. The grid size selected by the author is Δx=1320 m=1.32 km along the longitude and Δy=2600 m=2.6 km along the latitude. The uncertainty of earthquake location is about 2 km inland and about 3 km offshore. The values of Δx and Δy are both shorter than the location uncertainty. Can the simulated results based on the two values be applied to interpret the temporal variation and spatial distribution of aftershocks of the 2018 Hualien earthquake?
(4)For the 2D rectangle, its western part is on land, while its eastern part is offshore. The western part is below high mountains, while the eastern part is underwater. Hence, the vertical loading on the 2D rectangle is higher in the western part than in the eastern part. Hsu et al. (2025) showed that in southwestern Taiwan the excess seismicity rate is positively correlated with reduced NW-SE compression and/or decreasing vertical loading. This indicates the influence caused by vertical loading on seismicity. In the manuscript, the author wrote: ‘The Hualien M_w6.4 mainshock occurred in the complex convergence boundary of the Philippine Sea Plate and Eurasian Plate.’ This means that the geological structures and the values of physical parameters should be different in the two tectonic regimes. This would produce the difference in the stress diffusion between the two parts. Could the author clearly describe such differences due to different vertical loading.
[Reference]
Hsu, Y.J., R. Bürgmann, Z. Jiang, C.H. Tang, C.W. Johson, D.Y. Chen, H.H. Huang, M. Tang, and X. Yang (2025). Hydrologically-induced crustal stress changes and their association with seismicity rates in Taiwan, Earth Planet. Sci. Letts., 651, https://doi.org/ 10.1016/j.epsl.2024.119181.
(5)The mainshock showed thrust faulting and fault-to-fault jumping ruptures (Lee et al., 2019). Can the focal mechanism and rupture processes of the mainshock influence the spatial distribution of stresses and the evolution of stress diffusion?
[Reference]
Lee, S.J., T.C. Lin, T.Y. Liu, and T.P. Wong (2019). Fault‐to‐Fault jumping rupture of the 2018 M_w 6.4 Hualien earthquake in Eastern Taiwan. Seism. Res. Letts., 90(1). 30-39. https://doi.org/10.1785/0220180182
(6)Figure 1 displays that the number of aftershocks was larger in the eastern part (offshore) than in the western part (inland). This might indicate that the number of faults, sub-surface geological structures, and the values of physical parameters are different in the two parts. However, it does not seem able to apply the spatial distributions of stress changes at three depths (i.e., Figures 5-7) to explain the occurrence of aftershocks.
(7)As mentioned by the author, there are three faults in the study area. In fact, the number of faults is larger than three. Kuo-Chen et al. (2019) addressed a remarkable correlation between spatial distribution of aftershocks and the main fault along which the mainshock ruptured. Do the existence of those faults, particularly the main fault, which is not a straight line on the ground surface (cf. Lee et al., 2019; Kuo-Chen et al., 2019), influence the evolution of stress diffusion?
[Reference] Kuo-Chen H., Z.K. Guan, W.F. Sun, P.Y. Jhong, and D. Brown (2019). Aftershock sequence of the 2018 M_w 6.4 Hualien earthquake in Eastern Taiwan from a dense seismic array data set. Seism. Res. Letts., 90 (1), 60-67. https://doi.org/10.1785/ 0220180233
(8)There are remarkable differences in sub-surface geological structures inferred from 3D tomography between the western part and eastern one (e.g., Rau and Wu, 1995; Ma et al., 1996; Kim et al., 2005; Wu et al., 2007; Kuo-Chen et al., 2012). Can such differences influence the evolution of Coulomb stress?
[References]
Rau, R.J. and F.T. Wu (1995). Tomographic imaging of lithospheric structures under Taiwan. Earth Planet. Sci. Lett., 133, 517-532, Doi:10.1016/0012-821X(95)00076-O.
Ma K.F., J.H. Wang, and D. Zhao (1996). Three-dimensional seismic velocity structure of the crustal and uppermost mantle beneath Taiwan. J. Phys. Earth, 44, 85-105.
Wu, Y.M., C.H. Chang, L. Zhao, J.B.H. Shyu, Y.G. Chen, K. Shieh, and J.-P. Avouac (2007). Seismic tomography of Taiwan: Improved constraints from a dense network of strong‑motion stations. J. Geophys. Res., 112, B08312, doi:10.1029/2007JB004983.
Kim, K.H., Chiu, J.M., Pujo, J., K.C. Chen, B.S. Huang, Y.H. Yeh, P. Shen (2005). Three-dimensional V_p and V_s structural models associated with the active subduction and collision tectonics in the Taiwan region. Geophys. J. Intern., 162, 204-220.
(9)In Section 3, the author described the reasons how to select the values of model parameters for numerical simulations of the evolution of Coulomb stress after the 2018 Hualien earthquake in Taiwan. The main reason is for preventing the model from converging. Although this is an acceptable reason, the author should explain whether the values of model parameters can meet regional geological and seismological structures.

Minor Problems
Compared with the above-mentioned problems, the followings are minor.
Problems with Figures:
(1)The quality of figures is not good enough for readers. For example, Figures 5-7 are not good for readers.
Problems with References:
(1)The format of cited references should follow the Journal’s rules.
(2)the author should cite two important articles shown below for regional tectonics.
[References]
Tsai, Y.B., T.L. Teng, J.M. Chiu, and H.L. Liu (1977). Tectonic implications of the seismicity in the Taiwan region. Mem. Geol. Soc. China, 2, 13-41.
Wu, F.T. (1978). Recent tectonics of Taiwan. J. Phys. Earth, 2 (Suppl.), S265-S299.
(3)Few cited references, e.g., Zavyalov et al. (2024), are not listed in the Section of References.
(4)The KPP (or Fisher-KPP) equation is a very important nonlinear physical equation and has been widely applied in biology, ecology, and combustion theory. It is better to add cited references on Line 73 where the equation appeared at first in the text. The author may cite one of the following two articles or others:
Kolmogorov, A., I. Petrovskii, and N. Piskunov (1991). A study of the diffusion equation with increase in the amount of substance, and its application to a biological problem. In V.M. Tikhomirov, editor, Selected Works of A.N. Kolmogorov I, 248-270. Kluwer 1991, ISBN 90-277-2796-1. Translated by V.M. Volosov from Bull. Moscow Univ., Math. Mech. 1, 1-25.
El-Hachem M., S.W. McCue, W. Jin, Y. Du, and M.J. Simpson (2019). Revisiting the Fisher-Kolmogorov-Petrovsky-Piskunov equation to interpret the spreading–extinction dichotomy. Proc. R. Soc. A, 475:20190378. http://dx.doi.org/10.1098/rspa.2019.0378.
Problems with English Writing:
(1)The English writing is good. Nevertheless, there are some typo errors.
(2)The abstract is not concise.
(3)‘The Central Weather Bureau’ has been renamed ’the Central Weather Administration’ for a few years.

Citation: https://doi.org/10.5194/egusphere-2025-3348-CC1
- AC1: 'Reply on CC1', Boi-Yee Liao, 26 Dec 2025
  
  We sincerely thank Prof. Wang for his detailed and constructive comments on our manuscript.
  
  Given the scope and depth of the issues raised, we have carefully reviewed and consolidated our responses to ensure that each point is addressed rigorously and coherently.
  
  Several major comments concern fundamental aspects of model formulation, parameter interpretation, and physical–mathematical consistency, which will be fully and systematically addressed in the revised manuscript, along with detailed point-by-point responses.
  
  We appreciate the opportunity to improve the manuscript through this open discussion and thank Prof. Wang again for his valuable input.
  
  Citation: https://doi.org/10.5194/egusphere-2025-3348-AC1
RC1:
'Comment on egusphere-2025-3348', Anonymous Referee #1, 03 Jan 2026
Summary of Work:
In the present research, the author presented an iterative stress evolution model to the 2018 Mw 6.4 Hualien earthquake in Taiwan, focusing on depth-dependent Coulomb stress changes from 6 to 30 km. The model uses fixed parameters (α = 0.75, β = 0.7) and focuses the computation over on iterative counts for convergence, maximum stress value evolution, etc., and sensitivity analysis on parameters α (historical stress memory) and β (diffusion-reaction strength) using standard deviation (STD) and mean squared error (MSE) of stress fields at different depths.
The author beautifully connected the variation in STD to heterogeneity. For instance, STD increases initially with β, then decreases as diffusion smooths the field, which is linked to depths (6 km) show slower smoothing because of larger initial stress gradients, while deeper layers (30 km) preserve heterogeneity longer due to fewer iterations and weaker diffusion effects. Moreover, the persistence of heterogeneity is well-explained through α's role in retaining historical stress (e.g., higher α means slower decay, higher STD). Further, MSE complements this by measuring deviation from a "target" stable state—decreasing MSE with higher β shows how diffusion drives the system toward equilibrium, which is a solid proxy for persistence waning over time.
Moderate queries:
Question 1.
The structural features (faults) are missing in the figure 1. It can also include another parameter i.e. focal depth by filling the varying color in circle depicting earthquake location and magnitude.
Question 2.
Terms like n_critical, i.e., the iteration where diffusion dominates, are mentioned but not mathematically derived (Section 2.3). Further, a brief formula or sensitivity equation could strengthen why STD peaks differ by depth. These explanations feel descriptive rather than fully predictive and sufficient for intuition, which assists a better understanding without digging into earlier equations.
Question 3.
In my understanding, the STD and MSE are used effectively for model internals and the persistence is inferred from simulations, not directly compared to observed aftershock decay (e.g., Omori's law rates from the Hualien CatLog). Like aftershock patterns (e.g., clustering at 5–15 km), but section 3.2 could link STD trends more explicitly to empirical metrics like aftershock productivity or temporal decay exponents. I think, it would make the heterogeneity explanation more robust, especially since the model aims for aftershock prediction.
Question 4.
The STD captures spatial spread (heterogeneity), and MSE tracks convergence (persistence decay), but they do not fully address temporal aspects. For example, the paper notes iterations that are not directly related with time (shallow: hours/days per iteration; deep: months/years), i.e. it does not quantify time-scaling factors. This leaves persistence explanations a bit abstract—sufficient for parameter sensitivity, but could benefit from a time-calibrated example.
Question 5
STD and MSE trends with β at fixed α=0.75 (Figure 3a, b) and with α at fixed β=10 (Figure 3c, d) is explored thoroughly in result section. However, by exploring β variations across multiple α values, can we expect the following?
The peak STD occurs at different β (e.g., β=10 at 6 km vs. β=5 at 30 km with fix α=0.75). The higher α increases overall STD and potentially alters the optimum value of β at which diffusion dominates over reaction. Thus, under varying α, at shallower depths require higher β for smoothing, and could it highlight depth-specific couplings?

A large value of α slows stress decay, limiting the β's effective range to avoid divergence. Can we identify optimum parameter pairs by varying both parameters together that could optimize convergence rate without overshoot?

The α governs memory decay like R-S aging, while β scales diffusion/reaction like KPP terms (Table 1). Could the coupled variations better capture time-dependent behaviors (e.g., decreasing α(t) with increasing β(t) for aftershock decay) and provide a deeper understanding of stress redistribution in heterogeneous crust?

Is it possible that the extensions of analysis of explicitly 2D plot of β-α variations plot, may provide the refined model calibration, especially for real-time PSHA integration?

Question 6
The analysis assumes uniform parameter ranges (α=0.5–0.9, β=0.5–20) across depths, but real crust varies (e.g., rock rheology, fluids). The explanation is sufficient for the Hualien case, but generalizability to other tectonic settings (e.g., subduction zones) is not discussed, which might limit its perceived completeness.
The abstract mentions time-dependent α(t)/β(t), but the provided sections use fixed values—clarify if/when time-dependence is implemented, or emphasize fixed as an approximation.
Question 7
In section 3.3, the gradients are more predictive in the case of early (<50 days, <12 km) and absolute changes and backed by AUC/Molchan. A table summarizing these metrics across time/depth bins for: (i) static initial ΔCFF, (ii) dynamic ΔCFF(t), (iii) ∇σ(t), and (iv) combined. It may provide concrete evidence (e.g., AUC improvement of X% for gradients in early phase) and address why transitions occur.
Minor Query:
Banach Fixed-Point Theorem can be cited in line number 157.

The statement “Physically the contraction property reflects the Earth’s inherent damping mechanisms such as viscoelastic relaxation and afterslip” can be cited.

Please mention the figure number properly in line numbers 312 and 315.

Decision:
The statistical analysis strengthens the observations and potential claims for its operational use, i.e., PSHA. However, the article can be accepted after addressing the above moderate and minor queries raised.
Citation: https://doi.org/10.5194/egusphere-2025-3348-RC1
- AC3: 'Reply on RC1', Boi-Yee Liao, 08 Jan 2026
  
  The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2025/egusphere-2025-3348/egusphere-2025-3348-AC3-supplement.pdf
  
  Citation: https://doi.org/10.5194/egusphere-2025-3348-AC3
EC1:
'Comment on egusphere-2025-3348', Luciano Telesca, 07 Jan 2026

The feedback from CC1 and RC1 is both comprehensive and well-organized, offering clear insights into the challenges associated with the physical–mathematical modeling and interpretation. The authors are invited to revise their manuscript to fully address all the issues raised by the reviewers.

Citation: https://doi.org/10.5194/egusphere-2025-3348-EC1
- AC2: 'Reply on EC1', Boi-Yee Liao, 08 Jan 2026
  
  Dear Editor,
  Thank you for your message dated January 7 regarding the review of my manuscript (egusphere-2025-3348).
  I have now prepared the complete response documents, which thoroughly address all comments fromRC2.
  Best Regards
  Boi-Yee Liao
  
  Citation: https://doi.org/10.5194/egusphere-2025-3348-AC2

Peer review completion

AR – Author's response | RR – Referee report | ED – Editor decision | EF – Editorial file upload

AR by Boi-Yee Liao on behalf of the Authors (25 Jan 2026) Author's response Author's tracked changes Manuscript

ED: Referee Nomination & Report Request started (26 Jan 2026) by Luciano Telesca

RR by Anonymous Referee #1 (06 Feb 2026)

RR by Anonymous Referee #2 (18 Feb 2026)

ED: Publish as is (18 Feb 2026) by Luciano Telesca

AR by Boi-Yee Liao on behalf of the Authors (28 Feb 2026) Author's response Manuscript

Short summary

After major earthquakes, smaller shocks often follow, yet predicting where they will occur remains difficult. This study introduces a new method for tracking changes in underground stress after a large earthquake. Using the 2018 Hualien earthquake in Taiwan as a case study, we found that areas with strong stress differences provide clearer signals of future aftershocks than stress magnitude alone. This approach can improve short-term earthquake risk assessment and disaster response planning.

Beyond static forecasts: a dynamic stress gradient framework for high-resolution aftershock prediction and mitigation

Download

Interactive discussion

Peer review completion

Suggestions for revision or reasons for rejection