the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 20 Apr 2026
Hybrid implicit–explicit XFEM simulation of injection-induced seismicity: resolving multi-scale rupture nucleation and dynamics
Abstract. This study introduces a new hybrid implicit–explicit (IMEX) time-integration method for simulating injection-induced fault reactivation within a fully coupled hydromechanical extended finite element (XFEM) framework. Modeling these systems is inherently difficult because of the extreme difference in temporal scales – ranging from days for fluid diffusion to less than milliseconds for dynamic rupture. While fully implicit schemes provide unconditional stability, our results show that numerical accuracy during rupture nucleation is strictly controlled by the minimum time step used. Under-resolving these rapid transients causes delayed instability, lower slip rates, and underestimated seismic moments. Importantly, the temporal resolution needed to accurately capture rupture propagation is approximately the same as the critical time step set by the Courant–Friedrichs–Lewy (CFL) condition. To address this, the proposed IMEX strategy uses an implicit solver for quasi-static loading and switches smoothly to a fully dynamic explicit update once a slip-velocity threshold is reached. This hybrid approach significantly reduces computational costs – by 60–77 % compared to fully implicit simulations – while still accurately capturing fault slip, stress changes, and seismic magnitudes. Fully implicit solutions tend to predict slightly higher peak slip velocities and rupture speeds, but event timing and final deformation are consistent across both methods. Additionally, we show that the efficiency of the IMEX framework scales independently of the calculation of the complex slip tangent operator, providing a robust and efficient solution for multi-scale subsurface hazard assessment.
Competing interests: At least one of the (co-)authors is a member of the editorial board of Geoscientific Model Development.
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.- Preprint
(2473 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2026-1658', Anonymous Referee #1, 02 Jun 2026
-
AC1: 'Reply on RC1', mohammad Sabah, 12 Jun 2026
We sincerely thank Reviewer 1 for the careful evaluation of our manuscript and for the constructive and detailed comments. We appreciate the reviewer’s effort in identifying several parts of the manuscript where the presentation, definitions, mathematical notation, derivations, and explanations require further clarification. We agree that these comments will help improve the clarity, consistency, and readability of the manuscript. In the revised version, we will carefully address all comments raised by Reviewer 1 and update the manuscript accordingly.
Citation: https://doi.org/10.5194/egusphere-2026-1658-AC1
-
AC1: 'Reply on RC1', mohammad Sabah, 12 Jun 2026
-
RC2: 'Comment on egusphere-2026-1658', Anonymous Referee #2, 13 Jun 2026
Recommendation: Major revision
Overall assessment
The manuscript contains a potentially useful computational idea, but in its present form it substantially overstates the novelty, physical completeness, numerical performance, and practical applicability of the proposed method. I do not think that the current results are sufficient to support the main claims of the paper.
The formulation requires significant clarification, especially regarding the meaning of the proposed “IMEX” scheme, the governing poroelastic equations, the treatment of the fault, and the switch from slow loading to dynamic rupture. The numerical evidence is also not strong enough. In particular, the manuscript lacks independent validation of the new algorithm, systematic spatial and temporal convergence analysis, convincing demonstration of elastodynamic wave propagation, and performance tests at resolutions that would justify claims about scalability or operational applicability.
I therefore recommend major revision. However, I consider the required revision to be substantial rather than incremental. To become acceptable, the manuscript would need a much clearer and more defensible formulation, corrected or justified terminology, independent benchmark validation, convergence studies, larger and more realistic performance tests, and a more honest discussion of the limitations of the model. At present, the claims are too broad and too strong relative to what is actually demonstrated.
My main concerns are listed below.
The quasi-static to dynamic switching idea is not new by itself
The manuscript presents the transition from a quasi-static or implicit loading stage to a dynamic stage as a central novelty. I do not think this claim is sufficiently justified. The general idea of switching from slow quasi-static deformation to fully dynamic evolution once instability develops has already been used in related earthquake-cycle and rupture simulations. For example, Pampillón et al., “Dynamic and Quasi-Dynamic Modeling of Injection-Induced Earthquakes in Poroelastic Media”, already modeled injection-induced earthquake nucleation and dynamic rupture in poroelastic media, and Alkhimenkov et al. (2025), “Stress drop sequences in the simplest pressure-sensitive ideal elasto-plastic media: implications for earthquake cycles”, used a transition from slowly accumulated deformation to dynamic rupture evolution in an elasto-plastic setting.
1 The use of the term IMEX is not sufficiently justified
The manuscript presents the method as a hybrid implicit-explicit (IMEX) scheme. However, what is described appears to be closer to a solver-switching strategy, in which the simulation is advanced implicitly during the slow loading stage and then switched to an explicit mechanical update once a prescribed slip-velocity threshold is exceeded. This is not necessarily the same as a classical IMEX time discretization in the numerical-analysis sense, where implicit and explicit operators are usually combined within a single time-integration framework.
This distinction is important because the manuscript makes methodological claims based on the IMEX terminology. The authors should clarify whether the proposed method is a genuine IMEX discretization or a hybrid implicit-to-explicit switching algorithm. Moreover, the hydraulic equation appears to remain implicitly advanced during the nominally explicit phase. Therefore, the method does not seem to be fully explicit for the coupled hydromechanical problem. This should be stated clearly and reflected in the terminology.
2 The switching criterion is empirical, not convincingly physics-based
The switch from the implicit to the explicit phase is controlled by a prescribed slip-velocity threshold. I do not find this criterion sufficiently justified as a physics-based condition. A slip-velocity threshold can be a useful numerical trigger, but in the manuscript it is not derived from nucleation theory, process-zone size, characteristic slip distance, stress drop, wave speed, or any nondimensional physical argument.
The authors should either provide a clear physical derivation or interpretation of the threshold, or describe it more modestly as an empirical numerical switching parameter. Calling this switch “physics-based” is too strong in the present form.
3 Implicit schemes with large time steps are not only problematic for dynamic wave propagation
The manuscript discusses the limitation of implicit schemes mainly in relation to dynamic rupture and wave propagation. However, the need for small time steps is broader than that. During fluid injection, even before fully dynamic rupture, the system may develop rapid localization, accelerating slip, or new shear bands. These processes require adequate temporal resolution even if the governing problem is treated quasi-statically or implicitly.
Therefore, the issue is not only that implicit schemes become inefficient for elastodynamic wave propagation. Large implicit time steps may also miss or delay important transient processes during nucleation and localization. The manuscript should discuss this limitation more carefully.
4 The governing equations do not appear to represent full dynamic Biot poroelasticity
The manuscript refers to the formulation as dynamic poroelastic or elastodynamic Biot-type modeling. However, the governing equations presented appear to be a simplified u-p formulation rather than the full dynamic Biot system. In particular, I do not see the standard two-field inertial structure with the full 2 by 2 dynamic mass matrix, added-mass coefficients, tortuosity, and relative fluid acceleration terms.
The equations include an inertial correction in the Darcy-type pressure equation, but this is not the same as the full dynamic Biot poroelastic formulation. If the authors intentionally use a simplified formulation, they should state the approximation regime explicitly and justify it. Otherwise, the formulation should not be described as general dynamic Biot poroelasticity.
5 The model is restricted to reactivation of a prescribed fault
The manuscript is framed broadly in terms of injection-induced seismicity. However, the numerical model appears to address only reactivation of a pre-existing prescribed fault. It does not include bulk plasticity, damage, or a fracture-propagation criterion capable of generating new fractures or new shear bands in the surrounding medium.
This is an important limitation. In injection-induced seismicity, existing faults can indeed be reactivated, but new fractures and off-fault deformation can also develop depending on the stress state, injection conditions, and material properties. Since the present model cannot capture these processes, the authors should clearly restrict their claims to prescribed-fault reactivation. Broader statements about general injection-induced seismicity and physics-based hazard assessment should be moderated.
6 Rate-and-state friction is phenomenological and should not be oversold as physics-based
The manuscript uses a rate-and-state friction law. This is a valid and widely used constitutive model, but it is phenomenological. It is based on fitted parameters and an assumed state evolution law. Therefore, the use of rate-and-state friction alone does not make the model fully physics-based.
The authors should be more precise in their language. They should also provide all required details of the frictional initialization, including the reference velocity, reference state, initial state variable, initial slip velocity, and any regularization used at very small slip rates. These details are essential for reproducibility and for interpreting the switch based on slip velocity.
7 The initial stress and fault state need clearer specification
The manuscript provides initial stress values, but the initial mechanical state of the fault is not described in sufficient detail. The authors should specify the initial normal and shear tractions on the fault, the initial Coulomb ratio, the initial state variable, the initial slip velocity, and whether the model is equilibrated before injection begins.
This is important because the onset of fault reactivation depends strongly on the initial stress state and on how close the fault is to failure. Without this information, it is difficult to assess whether the reported nucleation behavior is physical or partly a consequence of initialization.
8 The treatment of the fault as a lower-dimensional interface needs more justification
The fault is represented as a lower-dimensional interface with hydraulic aperture and cubic-law flow. This is a useful modeling idealization, but it introduces assumptions that must be discussed more carefully.
The authors should clarify the effective fault thickness, how flow perpendicular to the fault plane is treated, how fluid exchange between fault and matrix is modeled, and how aperture-dependent permeability is coupled to the governing equations. If permeability, aperture, or other fault properties evolve during deformation, these evolution laws should be stated explicitly and consistently.
Introducing deformation-dependent hydraulic properties also makes the governing system nonlinear and modifies the standard poroelastic equations. The manuscript should reference previous work using similar modifications and should discuss the consistency and limitations of the resulting coupled formulation.
9 Fault geometry is unclear
The schematic figure and the computational setup appear to suggest different fault geometries. In one figure the fault/aperture representation looks elliptical or lens-like, whereas in the computational setup the fault appears as a long rectangular or line-like feature. The authors should clearly distinguish between the conceptual aperture schematic and the actual fault geometry used in the simulations.
The exact fault geometry, length, orientation, aperture distribution, and representation in the XFEM discretization should be stated clearly.
10 Dynamic wave propagation is not sufficiently demonstrated
The manuscript emphasizes dynamic rupture and elastodynamics, but I do not see clear snapshots of wave propagation generated by fault nucleation. If inertia and elastodynamic wave propagation are central to the method, the authors should show the resulting wavefield, for example snapshots of velocity, displacement, stress, or pressure perturbations, or receiver seismograms at several locations.
Without such evidence, it is difficult to assess whether the dynamic elastodynamic component of the solution is actually resolved or whether the analysis is mainly based on fault slip variables.
11 The mesh is too small to support strong performance and scalability claims
The numerical model uses a 200 by 200 mesh, corresponding to only 40,000 elements. This is a modest 2D problem. I do not think such a calculation can support strong claims about computational scalability, operational traffic-light systems, or future field-scale forecasting.
For comparison, existing high-resolution 2D and 3D numerical studies in related coupled mechanical and hydromechanical problems use more than 100 million voxels. A 10000 by 10000 two-dimensional voxel grid already contains 100 million elements, and modern three-dimensional GPU-based studies can reach hundreds of millions of cells, for example 500^3 = 125 million voxels. Against this background, a 200^2 simulation is far too small to demonstrate scalability or high-performance capability.
The authors may report that their hybrid strategy reduces CPU time relative to their own fully implicit implementation for this specific test. That is useful. However, it should not be presented as evidence of general scalability or readiness for operational hazard systems. To support performance claims, the authors should provide problem sizes, number of degrees of freedom, hardware, linear solver, preconditioner, number of time steps, Newton iterations, assembly time, solve time, and scaling behavior with increasing resolution.
12 A systematic convergence study is missing
A major weakness of the manuscript is the absence of a systematic spatial and temporal convergence analysis. The authors should run the solver at several spatial and temporal resolutions, for example four or five levels, and demonstrate convergence of key quantities: rupture time, peak slip velocity, final slip, rupture speed, seismic moment, pressure evolution, and wavefield measures.
Without such a study, it is unclear whether differences between the implicit and hybrid simulations are caused by the time-integration strategy, by spatial under-resolution, by the fault/contact discretization, or by insufficient temporal resolution. In dynamic rupture and rate-and-state friction problems, resolving the rupture process zone and nucleation length is essential. A CFL condition alone is not a sufficient accuracy criterion.
13 Validation of the new method is insufficient
The manuscript relies heavily on comparison with the authors’ own fully implicit solver. This is not sufficient validation of a new hybrid implicit-explicit switching algorithm. The fully implicit solver may be a useful reference, but the new method should be validated independently.
At minimum, the authors should include benchmark tests against analytical or semi-analytical solutions where possible. For the poroelastic part, classical tests such as Mandel-type problems or other fluid-injection solutions in poroelastic media would be appropriate. For the dynamic rupture part, comparison with established rupture benchmarks or highly resolved reference solutions is needed. The present validation does not convincingly establish the correctness of the coupled dynamic solver.
14 The conclusions overstate what has been demonstrated
The conclusions discuss broader applicability to 3D fault networks, realistic induced-seismicity forecasting, and adaptive traffic-light systems. These statements are not supported by the numerical evidence presented. The manuscript demonstrates a hybrid switching strategy on a small, homogeneous, two-dimensional, prescribed-fault problem. It does not demonstrate scalability to realistic 3D problems, complex fault networks, heterogeneous media, new fracture creation, or operational forecasting.
The authors should substantially narrow the conclusions and clearly separate what is demonstrated from what is only a possible future extension.
Citation: https://doi.org/10.5194/egusphere-2026-1658-RC2 -
AC2: 'Reply on RC2', mohammad Sabah, 13 Jun 2026
We thank the reviewer for the careful reading of our manuscript and for providing detailed and constructive comments. We acknowledge the reviewer’s concerns regarding the clarity of the formulation, terminology, numerical validation, convergence analysis, and the scope of the claims made in the manuscript. These comments are valuable and will help us improve the quality, rigor, and clarity of the revised manuscript.
In the revised version, we will carefully address the reviewer’s points by clarifying the numerical formulation and terminology, better explaining the assumptions and limitations of the proposed approach, moderating claims where necessary, and improving the discussion of validation, convergence, dynamic behavior, and computational performance. We will also revise the manuscript to more clearly distinguish between what is demonstrated in the present work and what remains a possible future extension.
We appreciate the reviewer’s constructive feedback and will revise the manuscript accordingly.
Citation: https://doi.org/10.5194/egusphere-2026-1658-AC2
-
AC2: 'Reply on RC2', mohammad Sabah, 13 Jun 2026
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 244 | 86 | 24 | 354 | 25 | 20 |
- HTML: 244
- PDF: 86
- XML: 24
- Total: 354
- BibTeX: 25
- EndNote: 20
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Dear authors, please find my comments on the manuscript in the attached file.