# Combined Boundary Element and Finite Methods for Modeling Fluid-Induced Seismicity in Fault Networks within Low-Permeability Rocks

Pierre Romanet, Jean-Paul Ampuero, Frederic Cappa, Marco Maria Scuderi, Stéphanie Chaillat

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Pierre Romanet, Jean-Paul Ampuero, Frederic Cappa, Marco Maria Scuderi, Stéphanie Chaillat. Combined Boundary Element and Finite Methods for Modeling Fluid-Induced Seismicity in Fault Networks within Low-Permeability Rocks. Geophysical Journal International, 2025, 243 (3), pp.ggaf377. 10.1093/gji/ggaf377. hal-05383298

HAL Id: hal-05383298

https://hal.science/hal-05383298v1

Submitted on 9 Feb 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Combined Boundary Element and Finite Volume
Methods for Modeling Fluid-Induced Seismicity in
Fault Networks within Low-Permeability Rocks

Pierre ROMANET¹,², Jean-Paul AMPUERO², Frédéric CAPPA²,³, Marco
Maria SCUDERI¹, and Stéphanie CHAILLAT⁴

¹Department of Earth Sciences, La Sapienza University of Rome, Rome, Italy

²Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, IRD, Geoazur, Sophia-Antipolis, France, email: romanet@geoazur.unice.fr

³Institut Universitaire de France, Paris, France

⁴Laboratoire POEMS, CNRS-INRIA-ENSTA Paris, Institut Polytechnique de Paris

February 9, 2026

Abstract

To better understand the mechanics of injection-induced seismicity, we developed a two-dimensional numerical code to simulate both seismic and aseismic slip on non-planar faults and fault networks driven by fluid diffusion along permeable faults, in an impervious host rock. Our approach integrates a boundary element method to model fault slip governed by rate-and-state friction with a finite-volume method to simulate fluid diffusion along fault networks. We demonstrate the capabilities of the method with two illustrative examples: (1) fluid injection inducing slow slip on a primary rough, rate-strengthening fault, which subsequently triggers microseismicity on nearby secondary, smaller faults, and (2) fluid injection on a single fault in a network
of intersecting faults, leading to fluid diffusion and reactivation of slip throughout the network. This work highlights the importance of distinguishing between mechanical and hydrological processes in the analysis of induced seismicity, providing a powerful tool for improving our understanding of fault behavior in response to fluid injection, in particular when a network of geometrically complex faults is involved.

**Keywords:**

Dynamics and mechanics of faulting TECTONOPHYSICS

Boundary element method GEOPHYSICAL METHODS

Induced seismicity SEISMOLOGY

Earthquake dynamics SEISMOLOGY

Mechanics, theory, and modelling TECTONOPHYSICS
## 1 Introduction

Understanding human-induced seismicity associated with anthropogenic activities has become an increasingly critical challenge for society. This urgency arises from the growing demand for energy, including green energy solutions, and the need to mitigate carbon dioxide (CO2) emissions. Many technologies developed to address these needs involve the injection and extraction of pressurized fluids into/from the subsurface. This is the case of Enhanced Geothermal Systems (EGS), which extract geothermal energy from deep rocks with initially low permeability, of disposal of wastewater from oil and gas operations into deep formations, and of carbon storage into sealed underground reservoirs. While these technologies provide substantial benefits in energy production and environmental management, they also introduce risks, particularly the potential of inducing earthquakes (Evans, 1966; Yerkes and Castle, 1976; Ellsworth, 2013; Grigoli et al., 2017; Keranen and Weingarten, 2018). In recent years, the rate of induced seismicity has increased significantly in regions with active fluid injections, leading to several relatively large earthquakes. Examples include the 2017 Mw 5.5 earthquake near an enhanced geothermal site in Pohang, South Korea (Grigoli et al., 2018; Kim et al., 2018; Lee et al., 2019). This earthquake caused substantial damage and has been proposed to be triggered by fluid diffusion from injection activities (Kim et al., 2018; Woo et al., 2019; Ellsworth et al., 2019). Another notable example is the 2011 Mw 5.7 and 5.0 earthquakes near Prague, in Oklahoma, USA (Keranen et al., 2013; Sumy et al., 2014), in which a first $M_w$ 5.0 earthquake was likely triggered by wastewater injection, and a subsequent $M_w$ 5.7 event was likely triggered by Coulomb stress changes. These events demonstrate that fluid injection can destabilize faults and trigger seismic events.

Fluid injection, however, does not systematically trigger a seismic rupture directly (Cornet, 2012; Ellsworth, 2013; Wei et al., 2015; Lengliné et al., 2017), injection-induced
aseismic slip is a process identified long ago (Hamilton and Meehan, 1971; Scotti and Cor-
net, 1994). During reservoir stimulations and field experiments, observations have shown
that aseismic slip can propagate from the injection point (Bhattacharya and Viesca, 2019),
potentially triggering seismicity further away (Guglielmi et al., 2015; Wei et al., 2015;
Cappa et al., 2019; Eyre et al., 2019). For instance, in hydraulic fracturing operations
within low permeability formations containing faults, aseismic slip may occur on stable
sections of a fault, which gradually load more distant, unstable regions. Eventually, this
loading can cause dynamic rupture when the aseismic slip front reaches a fault segment
with frictional rate-weakening behavior. This concept suggests that aseismic slip can prop-
agate from the reservoir depth to the seismogenic region of a fault, potentially triggering
earthquakes (Eyre et al., 2019; Atkinson et al., 2020; Eyre et al., 2022). Other cases, such
as those at the Soultz-sous-Forêts (France), Brawley (USA) and Cooper Basin (Australia)
geothermal sites, further illustrate how fluid injection can induce aseismic slip, causing
shear stress perturbations that trigger seismicity (Danré et al., 2024; Wang and Dunham,
2022). These cases highlight the role of fluid diffusion along fault systems in initiating
aseismic slip first, and then triggering seismicity.

In-situ experiments involving decametric fault reactivation have been essential for un-
derstanding the relationship between injected fluids, dilatancy and fault slip under con-
trolled conditions (Guglielmi et al., 2015; Duboeuf et al., 2017; Cappa et al., 2019; De Bar-
ros et al., 2024; Cappa et al., 2022a). For instance, recent observations at the Mont-
Terri underground laboratory (Switzerland) have shown that fluid injection in a highly
velocity-strengthening fault (Orellana et al., 2018; Bigaroni et al., 2023) with very low ini-
tial permeability (k = 10⁻¹⁷m²) can induce significant dilatant aseismic slip, which is a key
component of the deformation process in an interconnected fault network (Guglielmi et al.,
2020; Cappa et al., 2022b). This slip is often associated with an increase in low-frequency
tremor signals, linked to fluid-driven slip propagation (De Barros et al., 2023). In such
88

89

90

91

92

fault zones, structural heterogeneity, such as variations in fracture density, mineralogy, and
hydromechanical properties, plays a crucial role in controlling fluid diffusion, stress accu-
mulation and aseismic slip. Areas with rate-strengthening friction may experience aseismic
slip, while stress perturbations from aseismic slip may generate seismicity in regions with
rate-weakening properties (Cappa et al., 2019; De Barros et al., 2023).

93

94

95

96

97

98

99

100

101

Furthermore, centimeter scale laboratory experiments and decametric scale field tests have revealed the complexity of fault slip behavior with increasing fluid pressure (Scuderi and Collettini, 2016; Passelègue et al., 2018; Cappa et al., 2019; De Barros et al., 2023). Fluid pressurization can induce a variety of slip modes (seismic versus aseismic) in association with fault dilation, depending on fault hydromechanical properties, stress conditions, and fluid flow dynamics. These findings emphasize the importance of understanding the relationship between fault rheology, fluid-induced deformation, and stress changes due to fluid diffusion. These complexities underscore the need for models that capture the diversity of fault slip behaviors observed in experimental settings.

102

103

104

105

106

107

108

109

110

111

A significant body of research based on numerical modeling has focused on simulating
fluid-induced seismic and aseismic slip on faults governed by rate-and-state friction. Studies
primarily examined isolated planar faults with simplified fluid-fault coupling through the
effective normal stress, either under the quasi-dynamic approximation (McClure and Horne,
2011; Almakari et al., 2019; Dublanchet, 2019; Im and Avouac, 2025) or by including fully-
dynamic ruptures (Larochelle et al., 2021; Tan and Lui, 2025). Subsequent work expanded
this understanding by incorporating poroelastic effects, either in simplified representations
(Noda, 2022) or in more comprehensive forms (Norbeck and Horne, 2015; Heimisson et al.,
2022). These studies demonstrated the critical influence of bulk poroelastic properties,
hydraulic diffusivity, and friction evolution on rupture stability during fluid injection.

112

113

The impact of fluid pressure on complex fault geometries, including rough faults (Mau-
rer et al., 2020) and fault networks (Dieterich et al., 2015), has also been investigated.
These studies often relied on simplified assumptions, such as neglecting fluid flow or hy-
dromechanical coupling. Experimental evidence has shown that permeability is highly
sensitive to shear slip and normal stress, undergoing continuous changes throughout the
earthquake cycle (Barton et al., 1985; Esaki et al., 1999; Lee and Cho, 2002; Im et al., 2018,
2019). In particular, slip on non-planar faults induces changes in normal stress (Romanet
et al., 2020, 2024), which in turn can change permeability. In light of these observations,
McClure and Horne (2011) developed a rate-and-state friction model incorporating perme-
ability changes during earthquakes. Their work was later generalised to multiple non-planar
faults (Norbeck et al., 2016; Norbeck and Horne, 2016). Similarly, Cappa et al. (2018, 2019)
introduced a slip-dependent permeability model based on the cubic law to study the role
of aseismic slip during fluid injection in a rate-and-state fault. However, their approaches
were limited to modeling a single seismic event, as they neglected the permeability recovery.
Katsman and Ben-Avraham (2025), used a visco-elasto-plastic model, where permeability
change was incorporated through the change of volumetric strain, but it did not include
time dependent healing mechanism. These efforts emphasize the importance of modeling
fault permeability enhancement and dilatancy, particularly in low-permeability host rocks
such as shale layers and geothermal reservoirs. More recent work (Zhu et al., 2020; Yang
and Dunham, 2023; Ozawa et al., 2024) has addressed these limitations by incorporating
both permeability changes and recovery mechanisms, enabling the modeling of perme-
ability evolution throughout multiple earthquake cycles. The time scale of permeability
recovery is still debated and ranges from tens of minutes to years depending on the scale:
laboratory experiment or field observations (Elkhoury et al., 2006; Faoro et al., 2012; Xue
et al., 2013; Lai et al., 2014; Candela et al., 2014; Im et al., 2018). At the field scale, it
commonly ranges from months to ~ 2 years and could be an important factor for EGS
(Manga et al., 2012).

Motivated by natural cases, more realistic and complex fault network geometries need
to be incorporated in fluid-induced seismicity modeling. Fault systems in nature can be
difficult to image and exhibit complex activation mechanisms upon fluid injection. These
factors make the evaluation, prediction, and control of injection-induced earthquakes ex-
tremely challenging. To investigate fluid injection in fault networks, Discrete Fault Network
(DFN) models have been widely used to predict which faults might be reactivated, typically
under the assumption that the surrounding rock matrix is impermeable and fluid migra-
tion is confined to fault surfaces (Hicks et al., 1996; Kohl and Mégel, 2007; Maxwell et al.,
2015; Karvounis and Wiemer, 2022; Ciardo and Lecampion, 2023). While DFN models
offer valuable insights into potential fault reactivation, they often rely on the assumption
of a constant friction coefficient and Mohr-Coulomb failure, which limits their ability to
accurately capture the dynamics of slip evolution. Im and Avouac (2024) advanced this
approach by introducing a three-dimensional DFN rate-and-state model. However, their
model employed an analytical solution for fluid diffusion, which restricts the representation
of strong hydro-mechanical coupling, including critical effects such as permeability changes
due to deformation.

Each of these modeling approaches has advanced our understanding of fluid-induced fault slip in different ways: (1) by coupling fault slip with bulk poroelastic deformation (Norbeck and Horne, 2015; Heimisson et al., 2022; Noda, 2022), (2) by incorporating dynamic permeability changes and friction laws (McClure and Horne, 2011; Cappa et al., 2018, 2019; Zhu et al., 2020; Ozawa et al., 2024; Dunham, 2024) and (3) by capturing the complexities of fracture networks and geometry (Hicks et al., 1996; Kohl and Mégel, 2007; Maurer et al., 2020; Ciardo and Lecampion, 2023; Im and Avouac, 2024). Inevitably, each method has limitations due to either computational complexity, simplified representations or assumptions about fault permeability and connectivity.

In this study, we build upon recent advances by developing an integrated, computa-
tionally efficient approach for modeling fluid-induced seismicity in complex fault network
geometries and fluid injection scenarios. Our method enables the simulation of both aseis-
mic slip events and dynamic earthquakes during fluid injection. Fluid diffusion is modeled
along fault planes, with migration between interconnected faults. In Section 2, we intro-
duce our methodology, which couples an H-matrix accelerated Boundary Element Method
(BEM) for fault slip with a finite volume method (FVM) for fluid pressure diffusion. The
approach not only captures fault reactivation but also tracks the evolution of fault slip,
distinguishing between seismic and aseismic behavior. This distinction is particularly sig-
nificant, given the growing evidence that aseismic slip can outpace fluid diffusion, and
eventually trigger larger earthquakes at substantial distances from the injection point. In
Section 3, we showcase the method’s capabilities through two simulation examples that il-
lustrate how fluid-faulting interaction can control seismic activity in discrete fault networks
(DFN). Our modeling results are consistent with observed behaviors of induced earthquake
sequences, suggesting the key role of aseismic slip migration as a process controlling fluid-
faulting interaction leading to seismicity. Using these examples, we discuss the implications
of our results for understanding the interplay between fluid flow, aseismic slip and seismic
activity in fault networks within low permeability formations. The method developed in
this work can help gain insights that will offer valuable guidance to mitigate the risks of
induced seismicity in energy-related subsurface operations.

2 Method

2.1 Modeling fault slip coupled with fluid diffusion

We present a 2D hydro-mechanical model that integrates slip on a network of faults, gov-
erned by rate-and-state friction, with fluid diffusion along the fault surfaces. This model
is capable of capturing fault interactions within geometrically complex fault networks, and
of simulating seismic and aseismic slip. The model may be extended to accommodate the
non-linear evolution of fault permeability with ongoing slip. The main governing equations are described in this subsection, and the numerical methods used to solve them in the subsequent subsections.

The quasi-dynamic equilibrium of shear stress on the faults leads to the following equation that couples the mechanical and hydraulic components of our model:

$$ \underbrace{\tau_{\text{el}}(\Delta u) + \tau_{\text{load}} - \frac{\mu}{2\beta}V}_{\text{Fault shear stress}} = \underbrace{(\sigma_{\text{el}}(\Delta u) + \sigma_{\text{load}} - p)}_{\text{Effective normal stress}} \times \underbrace{f(V, \theta)}_{\text{Rate-and-state friction}} . \quad (1) $$

On the left-hand side, $\tau_{\text{el}}(\Delta u)$ is the static elastic shear stress due to the slip $\Delta u$ on each fault, $\tau_{\text{load}}$ the external shear stress loading on the faults, $V$ the slip velocity, $\mu$ the shear modulus, and $\beta$ the S-wave velocity. We adopt the quasi-dynamic approximation (Rice, 1993; Cochard and Madariaga, 1996), in which the radiation damping term $\frac{\mu}{2\beta}V$ represents the instantaneous dynamic shear stress reduction due to slip. It partially accounts for the elasto-dynamic effects due to wave propagation. On the right-hand side, the effective normal stress on the faults involves the static elastic normal stress due to slip $\sigma_{\text{el}}(\Delta u)$, the external normal stress load $\sigma_{\text{load}}$, and the pore fluid pressure $p$.

In principle, faults can open due to increasing pore fluid pressure if the effective normal stress vanishes. However, it remains difficult to model fault opening with rate-and-state friction because the state evolution equation has to be modified to disable state evolution during opening, and the ODE solver used here is not designed to handle boundary conditions involving inequalities related to fault opening. Some previous studies have accounted for changes in hydraulic aperture (McClure and Horne, 2011; Cappa et al., 2019) without loss of contact. In this article, the fault is not allowed to open, and we ensure that by stopping the calculation if the effective normal stress becomes tensile.
The friction coefficient is governed by rate-and-state friction: 211

$$f(V, \theta) = f_0 + a \log \left( \frac{V}{V_0} \right) + b \log \left( \frac{\theta V_0}{D_c} \right). \qquad (2)$$

It depends on the slip velocity $V$ and on a state variable $\theta$, and involves four parameters: a direct effect coefficient $a$, an evolution effect coefficient $b$, and a reference friction coefficient $f_0$ at a reference slip velocity $V_0$ (Marone, 1998). The evolution of the state variable is governed by the so-called aging law (though the model can accommodate other evolution laws): 212 213 214 215 216

$$\frac{d\theta}{dt} = 1 - \frac{V\theta}{D_c}, \qquad (3)$$

where $D_c$ is a critical slip distance. 217

The space and time evolution of fluid pressure on the faults is governed by the following diffusion equation, along each fault and across fault crossings: 218 219

$$\frac{\partial p}{\partial t}(\mathbf{x}, t) = \vec{\nabla} \cdot \left[ \alpha(\mathbf{x}) \vec{\nabla} p(\mathbf{x}, t) \right] + \frac{q(t)}{\phi(\beta_f + \beta_\phi)} \delta(\mathbf{x} - \mathbf{x}_{\text{inj}}), \qquad (4)$$

where $\vec{\nabla} \cdot$ is the divergence operator, $\vec{\nabla} p$ the pressure gradient, $\alpha(\mathbf{x}) = \frac{k(\mathbf{x})}{\phi\eta(\beta_\phi + \beta_f)}$ the diffusion coefficient, $k(\mathbf{x})$ the fault zone permeability, which may depend on the position along the fault, $\phi$ the porosity of the fault zone, $\eta$ the dynamic viscosity of the fluid, $\beta_\phi = \frac{1}{\phi} \frac{\partial \phi}{\partial P}$ the compressibility of the fault zone pores, $\beta_f = \frac{1}{\rho} \frac{\partial \rho}{\partial P}$ the compressibility of the fluid ($\rho$ is the fluid density), $q(t)$ the linear injection rate, $\delta(\mathbf{x})$ is the Dirac delta function and $\mathbf{x}_{\text{inj}}$ the injection point. This equation assumes that fault width is constant and that the specific storage $\phi(\beta_f + \beta_\phi)$ is constant along the faults. At junctions between two faults, diffusion occurs in all four possible directions. 220 221 222 223 224 225 226 227
## 2.2 Accelerated boundary element method for elastic stress interactions

The elastic shear stress $\tau_{el}(\mathbf{x})$ and normal stress $\sigma_{el}(\mathbf{x})$ at a given position $\mathbf{x}$ along a fault system due to a spatially distributed slip $\Delta u(\mathbf{x})$ is calculated using boundary integral equations (Tada and Yamashita, 1997; Romanet et al., 2020, 2024). For a non-planar fault system and in-plane slip (mode II), they are given by the following integrals along the fault network:

$$ \tau_{el}(\mathbf{x}) = \frac{\mu}{2\pi(1-\nu)} \int_{\text{faults}} K_{\text{grad}}^{t}(\mathbf{x}, \mathbf{y}) \frac{\partial}{\partial y^t} \Delta u(\mathbf{y}) dl(\mathbf{y}) + \frac{\mu}{2\pi(1-\nu)} \int_{\text{faults}} K_{\text{curv}}^{t}(\mathbf{x}, \mathbf{y}) \kappa^{t}(\mathbf{y}) \Delta u(\mathbf{y}) dl(\mathbf{y}), \quad (5) $$

$$ \sigma_{el}(\mathbf{x}) = \frac{\mu}{2\pi(1-\nu)} \int_{\text{faults}} K_{\text{grad}}^{n}(\mathbf{x}, \mathbf{y}) \frac{\partial}{\partial y^t} \Delta u(\mathbf{y}) dl(\mathbf{y}) + \frac{\mu}{2\pi(1-\nu)} \int_{\text{faults}} K_{\text{curv}}^{n}(\mathbf{x}, \mathbf{y}) \kappa^{t}(\mathbf{y}) \Delta u(\mathbf{y}) dl(\mathbf{y}), \quad (6) $$

where $\mu$ is the shear modulus and $\nu$ the Poisson's ratio of the host rock. The derivative of the slip along a fault is given by:

$$ \frac{\partial}{\partial y^t} \Delta u(\mathbf{y}) = \left( n_2(\mathbf{y}) \frac{\partial}{\partial y_1} - n_1(\mathbf{y}) \frac{\partial}{\partial y_2} \right) \Delta u(\mathbf{y}), \quad (7) $$

where $\mathbf{n}(\mathbf{y}) = (n_1, n_2)$ is the local unit vector normal to the fault. The absolute value of the local curvature along the fault $|\kappa^t(\mathbf{y})|$ can be interpreted geometrically as the inverse of the radius of the tangential circle to the fault. The stress transfer kernels $K_{\text{grad}}^t(\mathbf{x}, \mathbf{y})$, $K_{\text{curv}}^t(\mathbf{x}, \mathbf{y})$, $K_{\text{grad}}^n(\mathbf{x}, \mathbf{y})$ and $K_{\text{curv}}^n(\mathbf{x}, \mathbf{y})$ depend only on the fault network geometry, and are given in appendix A (see also eqs. C7 and C8 in Romanet et al. (2020)).

Equations (5) and (6) are discretised by decomposing the fault network geometry into a grid of non-overlapping linear straight elements and by assuming constant slip within
each element (Romanet et al., 2020). A constraint on the fault grid generation is that fault crossings must lie at the edge of neighboring elements, as shown in Fig. 1. Stresses are evaluated at the center of each element. After discretization, stresses are computed from slip values through a matrix-vector product between the boundary element matrix and the vector of slip values at all elements. This operation can be computationally expensive because the matrix is dense. To speed-up the stress evaluation, a data-sparse representation of the matrix is obtained with the hierarchical matrix technique (Hackbusch, 1999; Börm et al., 2003; Chaillat et al., 2017), see Ohtani et al. (2011); Bradley (2014); Sato and Ando (2021); Ozawa et al. (2023) for previous applications in earthquake modeling studies.

We have also implemented the antiplane (mode III) mode of slip. This is done by setting $\sigma_{el} = 0$ (there is no normal stress change in mode III) and $K_{curv}^t = 0$, and by changing the kernel $K_{grad}^t$ to the antiplane kernel given by Romanet et al. (2024).

## 2.3 Finite volume method for fluid diffusion along faults

The classic finite volume method (FVM) decomposes the fluid diffusion domain into control volumes $\Gamma_i$. In our case, considering diffusion along faults, each “volume” $\Gamma_i$ is a linear sub-segment of a fault, with length $L_i$. The primary variables in FVM are the averages of pressure $p$ and injection rate $q$ over each control volume, defined by

$$ 
\begin{aligned}
P_i(t) &= \frac{1}{L_i} \int_{L_i} p(\mathbf{x}, t) dl(\mathbf{x}), \\
Q_i(t) &= \frac{1}{L_i} \int_{L_i} q(\mathbf{x}, t) dl(\mathbf{x}).
\end{aligned} 
\qquad (8) $$

These averages are governed by the following spatially-discretized equations (Patankar, 2018):

$$ \frac{\partial P_i}{\partial t}(t) = \alpha_{i,i+1} \frac{P_{i+1}(t) - P_i(t)}{\Delta x_{i,i+1}^2} - \alpha_{i-1,i} \frac{P_i(t) - P_{i-1}(t)}{\Delta x_{i-1,i}^2} + \frac{Q_i(t)}{\phi(\beta_f + \beta_\phi)}, \qquad (9) $$
where $\alpha_{i,i+1} = \frac{2\Delta x_{i,i+1}}{L_{i+1}/\alpha_{i+1}+L_i/\alpha_i}$ is the harmonic average of the diffusivity between two neighboring elements $i$ and $i+1$, and $\Delta x_{i,i+1} = \frac{L_i+L_{i+1}}{2}$ is the distance between the center of elements $i$ and $i+1$.

The temporal discretization of the governing equations can be done by either explicit or implicit time schemes. Using an explicit scheme, which is conditionally stable, is not efficient because of the small time step required to ensure numerical stability. To circumvent this limitation, we adopt an implicit time scheme, the backward Euler scheme, which is unconditionally stable:

$$ \frac{P_i^{k+1} - P_i^k}{\Delta t} = \alpha_{i,i+1} \frac{P_{i+1}^{k+1} - P_i^{k+1}}{\Delta x_{i,i+1}^2} - \alpha_{i-1,i} \frac{P_i^{k+1} - P_{i-1}^{k+1}}{\Delta x_{i-1,i}^2} + \frac{Q_i^{k+1}}{\phi(\beta_f + \beta_\phi)}, \quad (10) $$

where the superscript $k+1$ stands for time $t_{k+1} = t_k + \Delta t$, where $\Delta t$ is the time step. It can be rewritten as:

$$ \alpha_{i-1,i} \frac{\Delta t}{\Delta x_{i-1,i}^2} P_{i-1}^{k+1} + \left( 1 - \alpha_{i,i+1} \frac{\Delta t}{\Delta x_{i,i+1}^2} - \alpha_{i-1,i} \frac{\Delta t}{\Delta x_{i-1,i}^2} \right) P_i^{k+1} + \alpha_{i,i+1} \frac{\Delta t}{\Delta x_{i,i+1}^2} P_{i+1}^{k+1} \\ = P_i^k + \Delta t \frac{Q_i^{k+1}}{\phi(\beta_f + \beta_\phi)}. \quad (11) $$

The discrete equation can be written in the form of a matrix-vector multiplication:

$$ \mathbf{A}\mathbf{P}^{k+1} = F(\mathbf{P}^k, \mathbf{Q}^{k+1}) \quad (12) $$

where $\mathbf{P}^k$ is the vector of pressure values at time $t_k$, $\mathbf{Q}^{k+1}$ is the injection rate at the new time $t_{k+1}$, $F$ is a linear function and $\mathbf{A}$ a tri-diagonal matrix that depends on $\Delta t$ and constant parameters. This system of linear equations has to be solved for $\mathbf{P}^{k+1}$ to update the fluid pressure at the new time $t_{k+1}$. In this work, it is solved with a conjugate gradient method.
<figure><img src="image_3.png" /><figcaption>Figure 1: The numerical grid for both Boundary Element Method (BEM) and Finite Volume Method (FVM). The connection between two faults is done at the element edges.</figcaption></figure>

The classic implicit FVM described above is for a single fault. We use a modification for two crossing faults, as required to handle diffusion along a fault network (Karimi-Fard et al., 2004; Berre et al., 2019). At the junction between two faults, diffusion occurs in all four possible directions. We further require fault junctions to lie at the edge of neighboring control volumes, as shown in Fig. 1. Let us consider the junction between faults 1 and 2, lying at the edge between elements i and i + 1 on fault 1 and between elements j and j + 1 on fault 2. To account for the effect of diffusion along fault 2, the FVM scheme for element i of fault 1 is updated to

$$ \frac{P_i^{k+1} - P_i^k}{\Delta t} = \underbrace{\alpha_{i,i-1} \frac{P_{i-1}^{k+1} - P_i^{k+1}}{(\Delta x_{i-1,i})^2}}_{\text{Effect of the left element}} + \underbrace{\alpha_{i,i+1} \frac{P_{i+1}^{k+1} - P_i^{k+1}}{(\Delta x_{i,i+1})^2}}_{\text{Effect of the right element}} + \underbrace{\frac{Q_i^{k+1}}{\phi(\beta_\phi + \beta_f)}}_{\text{Source term}} \\ + \underbrace{\alpha_{i,j} \frac{P_j^{k+1} - P_i^{k+1}}{(\Delta x_{i,j})^2}}_{\text{Effect of fault 2, element j}} + \underbrace{\alpha_{i,j+1} \frac{P_{j+1}^{k+1} - P_i^{k+1}}{(\Delta x_{i,j+1})^2}}_{\text{Effect of fault 2, element j+1}}. \qquad (13) $$

The first line in eq. (13) is the same as in a single fault. The second line accounts for the coupling of fluid pressure between the two faults. The other three elements attached to the fault crossing are updated in a similar way. With this modification, the matrix A is
no longer tri-diagonal.

## 2.4 Adaptive-time solver for the hydro-mechanical model

We adopt the same fault network discretisation for the BEM and FVM. We couple the FVM diffusion solver and the mechanical Ordinary Differential Equation (ODE) solver (Runge-Kutta 45) as follows:

1. An intended time step $\Delta t_{\text{try}}$ is suggested by the adaptive ODE solver Runge-Kutta 45, based on error estimates at 4th and 5th orders.

2. The time step is adjusted to guarantee the convergence of the pressure solver as follows. An attempt is made to compute the pressure $P^{k+1}$ at time $t_{k+1} = t_k + \Delta t_{\text{try}}$ using the FVM (equation 12). If within 200 iterations of the conjugate gradient solver the $L_\infty$ norm of the residual of equation 12 is not smaller than a prescribed tolerance $\epsilon$, $\Delta t_{\text{try}}$ is halved. This process is iterated until convergence to the desired accuracy is achieved. The resulting pressure is not used in subsequent calculations; this step is intended only to ensure convergence of the conjugate gradient method. The tolerance is set to be the same as the Runge-Kutta tolerance (here, we set an absolute tolerance of $\epsilon = 10^{-6}$ Pa by trial and error).

3. The ODE solver updates the slip velocity $V^{k+1}$, the state variable $\theta^{k+1}$ and the pressure $P^{k+1}$. In particular, the FVM solver is directly called to calculate the pressure at the required intermediate times, $P^{k+f}$ where $0 \le f \le 1$.

This time-solver can naturally deal with variable diffusivity that may depend on slip or on the effective normal stress, as long as the permeability evolution law is written as an ODE. Although not shown in the following examples, this is a feature that we plan to include later on in this code.
## 2.5 Code verification

Analytical solutions for fluid injection problems with rate-and-state friction are not available, thus we benchmark the code by comparison with other codes. We considered a community benchmark for fluid-induced aseismic slip, the SCEC SEAS benchmark problem BP06 (Lambert et al., 2025), which involves fluid injection on a planar fault governed by rate-and-state friction. The benchmark setup and results are described in Appendix B. The results show very good agreement between our FVM-BEM method, a discontinuous Galerkin method (Uphoff et al., 2023) and a spectral boundary element method (Dublanchet, 2019).

We further verified the convergence of our implementation of FVM by comparing numerical results obtained with different grid sizes to the analytical solution of the pressure diffusion scenario of the SCEC SEAS benchmark problem BP06 (Lambert et al., 2025). The results are reported in Appendix B. They show linear convergence as expected for the backward Euler scheme adopted here.

We also verified the convergence of our FVM-BEM in problems involving two crossing faults, by conducting simulations with different grid resolutions. The problem setting and results are presented in Appendix C). They demonstrate that our FVM-BEM reaches convergence upon grid refinement.

# 3 Examples of slip behavior and seismicity induced by fluid injection in fault networks

To showcase the capabilities of this seismo-hydro-mechanical model, we simulate two conceptual examples of fluid injection, each representative of realistic field-scale scenarios inspired by fault activation experiments at the decametric scale and reservoir stimulation
processes (e.g. Guglielmi et al. (2015)). The first example examines aseismic slip induced
by fluid injection on a rate-strengthening fault, which subsequently triggers microseismic-
ity on nearby secondary faults. The second example explores fluid diffusion and induced
aseismic slip within an interconnected fault network. In both scenarios, we adopt a simpli-
fied representation of hydraulic properties, assuming constant fault permeability (i.e., no
deformation-permeability coupling) and an impermeable surrounding rock (i.e., neglecting
poroelastic effects). The injected fluid is water at constant temperature (20°C). In the
following simulations, we pay particular attention to accurately resolving the process zone
size, given approximately by ~ 0.7L<sub>b</sub> (Ampuero and Rubin, 2008; Perfettini and Ampuero,
2008), where L<sub>b</sub> = (μD<sub>c</sub>) / ((1-ν)bσ′) and σ′ is the effective normal stress on the fault. The L<sub>b</sub> length
scale is discretized using at least 16 elements at the beginning of the simulations.

3.1 Microseismicity triggered near a rate-strengthening fault stim-
ulated by fluid injection

**3.1.1 Simulation settings**

The objective of this example is to demonstrate how fluid injection can induce aseismic slip on a rate-strengthening fault, which subsequently triggers microseismic activity on adjacent secondary rate-weakening faults through elastic stress transfer. We simulate localized fluid injection on a 10 m-long permeable fault with diffusivity $\alpha = 10^{-3}$ m$^2$/s, surrounded by a network of smaller, secondary faults. Values of model parameters are given in Tab. 1.

The main fault is represented as a self-similar geometrically rough fault (Power and Tullis, 1991; Renard et al., 2006; Candela et al., 2011; Dunham et al., 2011), with deviations from planarity characterized by an amplitude-to-wavelength ratio of $5 \times 10^{-4}$, and a minimum wavelength of $\lambda_{min} = 20$ cm. Surrounding the main fault are 800 secondary faults with equal length, each 100 times shorter than the main fault, distributed on a regular grid
that spans an area of 15 m by 8 m. They are oriented with variable angles, randomly drawn from a uniform distribution between -30° and +30° relative to the main fault's orientation. The hydraulic and frictional properties of the faults, as well as the elastic properties of the surrounding rock, are assumed to be uniform and set to typical values consistent with crustal reservoir conditions (Table 1). We impose a no-flow condition at the end of the main fault, so that the injected fluid remains trapped within the fault. The main fault has rate-strengthening frictional properties ($a > b$) and thus can slip aseismically. The secondary faults are velocity-weakening ($a < b$) and thus can rupture dynamically when the Coulomb stress changes caused by slip on the main fault are positive and sufficiently large. To ensure computational efficiency, each secondary fault is modeled with a single element, allowing us to capture the temporal evolution of their moment release while omitting the spatial details of their rupture propagation.

The simulation is designed to mimic a faulted crustal reservoir in an intraplate tectonic regime, where the external stressing rate is supposed to be very slow or inexistent (e.g. Townend and Zoback (2000)). In the following simulations, we do not impose any tectonic loading (Larochelle et al., 2021). This approach is also justified by the fact that stress changes due to fluid injection typically occur much more rapidly than those caused by tectonic loading. This is evidenced by the sharp increase in seismicity rates in the central US following the onset of wastewater injection (Ellsworth, 2013; Magnani et al., 2017). In this context, the fault slip history in our model is governed solely by the stress perturbations resulting from fluid injection and diffusion.

The initial conditions are such that all faults are initially sliding at steady state (i.e. $V_0\theta_0/D_c = 1$). The system is then allowed to evolve for 10 years to minimize the influence of initial conditions on the subsequent fluid injection. When injection starts the faults are below steady-state. Following this period of model initialization, fluid is injected at a constant rate ($q = 1.3416 \times 10^{-6} \text{ m s}^{-1}$) on the main fault, at a point 1 m away from
its center. Injection lasts 7,000 seconds and reaches a peak pressure of 2 MPa. In the simulation, the evolution of fluid pressure, fault slip, slip velocity, stress and friction are calculated. Slip is defined to be seismic when slip velocity exceeds a threshold of 1 mm/s. The slip velocity does not diverge, owing to the stabilizing effect of the radiation damping term of the quasi-dynamic approximation (Rice, 1993).

### 3.1.2 Injection-induced slow slip evolution

Before fluid injection, the main fault slides very slowly, with a maximum slip velocity of the order of $5.4 \times 10^{-15}$ m/s (Fig. 2.a). The onset of fluid injection initiates a slow aseismic slip transient that nucleates at the injection point (Fig. 2.b). The subsequent evolution of injection-induced slow slip on the main fault (Figs. 2 and 3) follows distinct stages similar to previous studies.

For a velocity-strengthening rate-and-state fault, injection-induced slow slip has been shown theoretically to follow a two-stage development (Dublanchet, 2019). In the early stage of injection, when the state evolution and the elastic stress transfer are negligible, the logarithm of slip velocity is approximately proportional to the fluid pressure change (equation A.2 of Dublanchet (2019)). This behavior is observed in our simulation during the early phase of slow slip (Fig., 3.b). It arises from the direct effect in rate-and-state friction, and hence is very different from what can be observed in models that assume constant friction (Bhattacharya and Viesca, 2019; Viesca, 2021; Sáez et al., 2022) or slip-weakening friction (Garagash and Germanovich, 2012; Sáez and Lecampion, 2024). In this early stage, the extension of the slow-slip front is mainly controlled by fluid diffusion.

In a second stage, the region of aseismic slip expands and takes the form of a growing quasi-static crack (Figs. 2.c,d and 3.c,d). We define the slow slip front by the position where slip velocity reaches a threshold of $10^{-11}$ m/s, a value much smaller than the peak slip velocity. We observe a nearly linear growth of the slow-slip front with time up to
approximately 1500 s after injection onset (Fig. 3.a). At later times, when the crack-like behavior is well developed, a second quasi-linear growth phase emerges, although perturbed by fault end effects. Theoretical models of fluid injection at constant pressure on a single planar fault, governed by constant friction (Bhattacharya and Viesca, 2019; Viesca, 2021; Sáez et al., 2022) or by slip-weakening friction (Garagash and Germanovich, 2012; Sáez and Lecampion, 2024), predict a diffusive slow slip migration in which the slip front grows proportionally to $\sqrt{t}$. However, in the case of injection at constant flow rate on a 1D fault, as in our simulation, linear propagation of slow slip is predicted by models with constant friction (appendix D in Sáez et al. (2022)) and with velocity-strengthening rate-and-state friction (Dublanchet, 2019). This linear crack-like growth is also predicted under rate-and-state friction by the fracture mechanics theory of Garagash (2021). We note, however, that this is a feature of the low dimensionality of our model: on a 2D rate-and-state fault, the propagation of injection-induced slow slip under constant injection flow rate is predicted to be diffusive rather than linear (Danré et al., 2024).

Comparing the patterns of slow slip propagation and fluid pressure diffusion provides insights on their hydro-mechanical interactions. While the slow slip front can be defined without major ambiguity, the usual concept of pressure diffusion front implies a selected pressure contour, whose choice is largely arbitrary. Indeed, the diffusion front is often defined by

$$x_{\text{diff}} = \sqrt{\alpha t}, \qquad (14)$$

where $x_{\text{diff}}$ is the distance between the front and the injection point and $t$ is the elapsed time since the start of the injection (Shapiro et al., 2002; Bhattacharya and Viesca, 2019; Ciardo and Lecampion, 2023). For a point injection at constant pressure $P_0$ along a 1D fault, using the analytical solution for this 1D diffusion problem (Carslaw and Jaeger, 1959) we find that this front definition corresponds to the pressure contour $P(x, t) = 0.157 P_0$.
Some authors include a factor inside the square root in equation 14 (e.g., Shapiro et al. (1997)), which changes the pressure contour associated to the pressure front. For injection at constant flow rate, as in our simulations, the fluid pressure is (e.g. appendix D in Sáez et al. (2022)):

$$P(x, t) = P_{\max}(t) \Pi(\xi) \qquad (15)$$

where $\xi = x/\sqrt{4\alpha t}$ is a self-similar variable, and

$$P_{\max}(t) = \frac{q_0 \eta}{\sqrt{\pi k}} \sqrt{\alpha t} \qquad (16)$$
$$\Pi(\xi) = \exp(-\xi^2) - \sqrt{\pi} |\xi| \text{Erfc}(|\xi|)$$

For this case, Sáez et al. (2022) defined the diffusion front as

$$x_{\text{diff}} = \sqrt{4\alpha t}. \qquad (17)$$

This definition corresponds to $\xi = 1$ and hence to the following contour of the normalized fluid pressure: $P(x_{\text{diff}}, t)/P_{\max}(t) = \Pi(\xi = 1) = 0.089$. Clearly, a different choice of the factor inside the square root in equation 17 would associate the pressure front to a different normalized-pressure contour, revealing a degree of arbitrariness in the definition of the pressure front (Sáez et al., 2025). As illustrated in Fig. 3.a, depending on the chosen pressure contour, different interpretations arise regarding whether slow slip outpaces pressure diffusion. If the pressure front is defined by the 0.089 $P_{\max}(t)$ contour (bold red line), consistent with Sáez et al. (2022), the slow slip lags behind the pressure front until ~2000 s, then outpaces it. If the pressure front is defined by a higher pressure contour, such as 0.2 $P_{\max}(t)$, slow slip outpaces it at an earlier time, after around 1500 seconds after injection begins. In contrast, if we define the pressure front by the 0.01 $P_{\max}(t)$ contour, slow slip does not outpace pressure diffusion. This underscores the need for
careful interpretation when comparing the spatial evolution of fluid pressure and fault slip.

### 3.1.3 Microseismicity triggered by slow slip

The aseismic slip on the main fault transfers stresses to the nearby secondary faults, causing their reactivation and triggering micro-seismicity (Figure 2.e,f). In our simulations, seismicity does not start right after the onset of injection. Instead, two clusters of seismicity, one on each side of the main fault, start only after the slow slip reaches the corresponding tip of the fault, which happens before injection ceases (the left tip is reached ~ 45 minutes after injection starts, Figure 3.a; see also the video described in the Data Availability Statement). As fluid pressure continues to diffuse within the main fault, the evolving slow-slip causes increased Coulomb stress at the fault tips, which initiates microseismicity on the adjacent, rate-weakening secondary faults. Seismicity remains active even after injection has stopped. As diffusion proceeds, a constant fluid pressure is achieved along the fault (Figs. 2.f). Continued aseismic slip on the primary fault continues to load the surrounding secondary faults, thereby sustaining seismicity, which continues to propagate outward from the main fault tips (Fig. 2.f).

In this model, there is no fluid diffusion from the main fault to the secondary faults nor poroelastic effects that could contribute to their failure. As expected, the location of seismicity corresponds to the regions of positive Coulomb stress changes concentrated at the edges of the main fault. There is stress feedback between the main fault and the secondary faults, but this effect is too small to change the dynamics of the main fault, as shown in appendix D by comparison to a simulation without secondary faults.

To better understand the evolution of microseismicity, and its relation with fluid pressure diffusion and aseismic slip, we analyzed its time-dependent migration. This is particularly important because field observations have shown that injection-induced aseismic deformation can alter local stress fields, controlling the migration of seismicity and, in some
cases, reactivating misoriented faults far from the injection point (Schoenball et al., 2014; Duboeuf et al., 2021; Goebel et al., 2017). To confirm the robustness of the observed seismicity patterns, we performed additional simulations in which we randomly changed the orientation of the secondary faults and the roughness of the main fault (see also Appendix D), while maintaining all other parameters unchanged. During injection, the spatiotemporal distribution of seismicity in all these simulations (Figure 4.a) exhibits an approximately linear migration away from the fault. This is in contrast with a diffusion-like pattern.

The post-injection seismicity also has interesting relations to the slow slip. The event most distant from the main fault occurs shortly after the end of injection. Afterwards, no further migration of seismicity is observed, and the seismicity rate decreases rapidly, although it does not vanish over the duration of our simulations (Fig. 4.b). Thus, slow slip induced by fluid injection perturbs the secondary faults persistently so that, after injection stops, seismicity continues for a long time compared to the injection duration. Indeed, after injection ends, the seismicity rate decays following approximately a $1/t$ Omori law and is roughly proportional to the mean slip velocity on the main fault (inset of Fig. 4.b). This behavior is analogous to the aftershock model of Perfettini and Avouac (2007), in which a velocity-strengthening fault responds to a co-seismic static stress perturbation by a long-term Omori-type slip velocity decay (e.g. Perfettini and Ampuero, 2008), which in turn triggers off-fault seismicity. The proportionality between the seismicity rate and the loading rate imposed by slip velocity is what would be obtained in a simple Coulomb failure model, but is also expected in a population of rate-and-state faults at time scales longer than their characteristic nucleation time (Dieterich, 1994).
<table><caption>Table 1: Parameter values for the simulation of micro-seismicity triggered by fluid injection.</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td colspan="3">Frictional properties of the main fault</td></tr><tr><td>a</td><td>0.006</td><td>Direct effect</td></tr><tr><td>b</td><td>0.005</td><td>Evolution effect</td></tr><tr><td>D<sub>c</sub></td><td>1 µm</td><td>Critical Distance</td></tr><tr><td>f<sub>0</sub></td><td>0.6</td><td>Reference friction</td></tr><tr><td>V<sub>0</sub></td><td>1 × 10<sup>-6</sup> m s<sup>-1</sup></td><td>Reference velocity</td></tr><tr><td colspan="3">Frictional properties of the secondary faults</td></tr><tr><td>a</td><td>0.005</td><td>Direct effect</td></tr><tr><td>b</td><td>0.007</td><td>Evolution effect</td></tr><tr><td>D<sub>c</sub></td><td>4 × 10<sup>-3</sup> µm</td><td>Critical Distance</td></tr><tr><td>f<sub>0</sub></td><td>0.6</td><td>Reference friction</td></tr><tr><td>V<sub>0</sub></td><td>1 × 10<sup>-6</sup> m s<sup>-1</sup></td><td>Reference velocity</td></tr><tr><td colspan="3">Rock elastic properties</td></tr><tr><td>μ</td><td>3.2040 × 10<sup>10</sup> Pa</td><td>Shear modulus</td></tr><tr><td>ν</td><td>0.1541</td><td>Poisson's ratio</td></tr><tr><td colspan="3">Fault hydraulic and fluid properties</td></tr><tr><td>k</td><td>1 × 10<sup>-15</sup> m<sup>2</sup></td><td>Permeability</td></tr><tr><td>β<sub>f</sub></td><td>9 × 10<sup>-9</sup> Pa<sup>-1</sup></td><td>Compressibility of the fluid</td></tr><tr><td>β<sub>φ</sub></td><td>1 × 10<sup>-9</sup> Pa<sup>-1</sup></td><td>Compressibility of the pore</td></tr><tr><td>φ</td><td>0.1</td><td>Porosity</td></tr><tr><td>η</td><td>1 × 10<sup>-3</sup> Pa s</td><td>Dynamic viscosity</td></tr><tr><td>q</td><td>1.3416 × 10<sup>-6</sup> m s<sup>-1</sup></td><td>Flow rate of injection</td></tr><tr><td>t<sub>inj</sub></td><td>10 years</td><td>Beginning of injection</td></tr><tr><td>Δt<sub>inj</sub></td><td>7000 s</td><td>Duration of injection</td></tr><tr><td colspan="3">Initial values for mechanical conditions</td></tr><tr><td>V<sub>ini</sub></td><td>10<sup>-12</sup> m s<sup>-1</sup></td><td>Initial slip velocity</td></tr><tr><td>θ<sub>ini</sub></td><td>10<sup>6</sup> s</td><td>Initial state</td></tr><tr><td>σ<sub>load</sub></td><td>5 MPa</td><td>Normal stress</td></tr><tr><td>τ<sub>ini</sub></td><td>2.59 MPa</td><td>Initial shear stress on main fault</td></tr><tr><td>τ<sub>ini</sub></td><td>2.65 MPa</td><td>Initial shear stress on secondary faults</td></tr><tr><td colspan="3">Length scales</td></tr><tr><td>ds</td><td>5 mm</td><td>Grid size</td></tr><tr><td>L<sub>b</sub></td><td>1.3329 m</td><td>Process zone size on main fault (Ampuero and Rubin, 2008)</td></tr></tbody></table>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 2: Evolution of injection-induced slow slip on a rough main fault and triggered earthquakes on secondary faults. In each panel, the geometry of the fault network (comprising one main fault and multiple secondary faults) is shown at the bottom, the spatial distribution of fluid pressure is shown along the main fault (color scale on the bottom right), slip velocity profiles are shown as curves above the main fault (in logarithmic scale, vertical axis), secondary faults that have slept seismically are indicated by red stars, and the 2D spatial distribution of Coulomb stress along slip planes parallel to the main fault is shown in colors at the bottom (color scale on the bottom left). The selected snapshot times for each panel, relative to the onset of injection, are: (a) -25 seconds (right before injection), (b) 5 minutes, (c) 20 minutes, (d) 30 minutes, (e) 2 hours and (f) 1 day.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 3: (a) Slow-slip front (blue curve) and selected contours of normalized fluid pressure (normalized by the maximum pressure at each time, red curves; the darker red curve corresponds to the pressure front defined consistently with equation 17). (b-d) Spatial distribution of slip velocity (blue) and fluid pressure (red) at three selected times indicated in (a) by vertical dashed lines: (b) when the logarithm of the slip velocity is approximately proportional to the fluid pressure, (c) when the elastic effects start to dominate and slip velocity is still localized, and (d) when the slow slip shows crack-like expansion.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 4: (a) Space-time distribution of seismicity. Time is relative to the start of injection and distance is relative to the closest point on the main fault. Stars show the micro-seismicity of (red) the simulation shown in Fig. 2 and (light yellow) five additional simulations with different random realizations of fault roughness and secondary fault orientations. (b) Seismicity rate of the composite catalog of the 6 simulations, calculated in bins of ±15 minutes, and mean slip velocity of the main fault from the simulation of Fig. 2. The inset shows a zoomed-in view of the same quantities, re-scaled to highlight the approximate proportionality between the seismicity rate and the slip-velocity after the end of injection.</figcaption></figure>

## 3.2 Fluid diffusion and aseismic slip propagation in a fault network

In the second example, we examine the response of a hydraulically connected network of faults to fluid injection. The geometry of the fault network is defined in Fig. 5 and parameter values are listed in Tab. 2. Unlike the first example, all the faults in this scenario are planar, have comparable lengths, and the rupture on each fault is spatially well-resolved with grid sizes smaller than the process zone size (Rubin and Ampuero, 2005). The frictional properties are uniform across all faults, exhibiting a rate-strengthening behavior.
Slip occurs in mode II (in plane shear), with interactions between faults facilitated by static stress transfer. The injection procedure involves constant injection rate at 2 × 10⁻⁶ m/s over a 2000 seconds time interval, pressurizing the fluid up to 1.6 MPa. As injection begins, fluid diffusion occurs throughout the fault network. There are no-flow conditions at the end of the faults, hence all the injected fluid remains in the fault network.

<figure><img src="image_12.png" /><figcaption>Figure 5: Geometry of the fault network. The angle between the two families of faults is 60°.</figcaption></figure>

The response of the fault network to fluid injection is visualized through the snapshots of slip velocity and fluid pressure shown in 2D view in Fig. 6 and as a function of distance from the injection point in Fig. 7. Before injection, all faults are slowly slipping due to a small applied shear stressing rate of 0.001 Pa/s (Figs. 6.a and 7.a). Fluid injection is initiated on Fault 10. Shortly after, fluid pressure at the injection point increases rapidly, triggering a slow slip event along the injection fault within 5 minutes (Fig. 6.b). In these early stages, there is an approximately linear relation between the rise in fluid pressure and
the logarithm of slip velocity on Fault 10 (Fig. 7.b), as expected in rate-and-state models when the changes in slip and state variable are still small (Dublanchet, 2019; Larochelle et al., 2021). As the injection ends, 30 minutes after it started, a complex slip history develops across the fault network (Fig. 6.c and 7.c). Notably, while fluid diffusion remains largely confined within the injection fault (Fault 10), the slow slip propagates further, reaching Faults 9 and 5 (Figs. 6.c). From Fig. 7.c we can discriminate that the slow slip on Fault 5 is driven by the redistribution of stress across faults rather than by the fluid diffusion process. Indeed, fluid pressure suddenly drops well before reaching Fault 5, at the intersection between Fault 10 and Fault 9, as fluid diffusion transitions from a single pathway to three directions. Well after injection has ended, after 12 hours, the fluid pressure has decreased to a minimum value and is slowly diffusing along the fault network (Fig. 6.d and 7.d). Nonetheless, the fault network is still quite active: faults far from the injection point (i.e., Faults 5, 7 and 8) are still slowly slipping while the injection fault (Fault 10) and the first connected fault (Fault 9) are decelerating.

To better characterize the slow slip propagation, we define the position of the slow slip front by a slip-velocity threshold of $10^{-11}$ m/s. We plot its distance from the injection point in Fig. 8, together with different contours of normalized fluid pressure $P(x,t)/P_{max}(t)$ (red lines). The acceleration of the normalized-pressure contour after the end of injection is an artifact of the decreasing $P_{max}(t)$. As in the previous example, different choices of a contour lead to different interpretations of when slow-slip outpaces the diffusion front. The slow slip front is ahead of the nominal contour of $0.089P_{max}(t)$ from around 600 s after the injection.

Something noticeable is that the slow-slip can propagate more easily through the fault network than the fluid pressure. This is because at each fault crossing, the fluid pressure has to diffuse to four directions, which impedes its propagation. This is a major difference from the previous simulation, in which the fluid was diffusing only along a single fault.
547

The apparent cessation of fluid pressure diffusion observed in figure 8, approximately
two meters from the injection point, is attributed to the closed end of fault 10. In fact, this
closed end is filled with fluid and the only remaining path for fluid to diffuse is through
the fault network. This observation suggests that fluid diffusion within a closed network of
faults, under the simplified scenario of constant diffusivity, exhibits an apparent diffusivity
that depends on the injection location, the dead-ends of fault and the connectivity between
faults.

548

549

550

551

552

553

554

This example highlights the complexity of the interaction between fluid diffusion and
stress transfer due to slow slip transients. It is important to note that this is a very
simplified case scenario and that more sophisticated 3D geometry, heterogeneous frictional
properties and injection history may give rise to the level of complexity that is observed
in nature.

555

556

557

558
<table><caption>Table 2: Parameter values for the simulation of fluid diffusion and aseismic slip propagation along a fault network.</caption><thead><tr><th colspan="3">Faul frictional properties</th></tr></thead><tbody><tr><td>a</td><td>0.005</td><td>Direct effect</td></tr><tr><td>b</td><td>0.004</td><td>Evolution effect</td></tr><tr><td>D<sub>c</sub></td><td>1 µm</td><td>Critical Distance</td></tr><tr><td>f<sub>0</sub></td><td>0.6</td><td>Reference friction</td></tr><tr><td>V<sub>0</sub></td><td>1 µm s<sup>-1</sup></td><td>Reference velocity</td></tr><tr><th colspan="3">Elastic properties of rock</th></tr><tr><td>μ</td><td>32.040 GPa</td><td>Shear modulus</td></tr><tr><td>ν</td><td>0.1541</td><td>Poisson's ratio</td></tr><tr><th colspan="3">Fault hydraulic and fluid properties</th></tr><tr><td>k</td><td>1 × 10<sup>-15</sup> m<sup>2</sup></td><td>Permeability</td></tr><tr><td>β<sub>f</sub></td><td>9 × 10<sup>-9</sup> Pa<sup>-1</sup></td><td>Compressibility of the fluid</td></tr><tr><td>β<sub>φ</sub></td><td>1 × 10<sup>-9</sup> Pa<sup>-1</sup></td><td>Compressibility of the pore</td></tr><tr><td>φ</td><td>0.1</td><td>Porosity</td></tr><tr><td>η</td><td>1 × 10<sup>-3</sup> Pa.s</td><td>Dynamic viscosity</td></tr><tr><td>q</td><td>2 × 10<sup>-6</sup> m s<sup>-1</sup></td><td>Flow rate of injection</td></tr><tr><td>t<sub>inj</sub></td><td>10 years</td><td>Beginning of injection</td></tr><tr><td>Δt<sub>inj</sub></td><td>2000 s</td><td>Duration of injection</td></tr><tr><th colspan="3">Initial values for mechanical conditions</th></tr><tr><td>V<sub>ini</sub></td><td>1 × 10<sup>-12</sup> m s<sup>-1</sup></td><td>Initial slip velocity</td></tr><tr><td>θ<sub>ini</sub></td><td>10<sup>6</sup> s</td><td>Initial state</td></tr><tr><td>σ<sub>load</sub></td><td>5 MPa</td><td>Normal stress</td></tr><tr><td>τ<sub>ini</sub></td><td>2.65 MPa</td><td>Shear traction</td></tr><tr><th colspan="3">Length scales</th></tr><tr><td>ds</td><td>10 cm</td><td>Grid size</td></tr><tr><td>L<sub>b</sub></td><td>1.67 m</td><td>Process zone size (Ampuero and Rubin, 2008; Perfettini and Ampuero, 2009)</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Figure 6: Snapshots of slip velocity and fluid pressure along the fault network. Selected times, relative to the onset of injection, are: (a) −4 minutes 20s (right before injection starts), (b) 5 minutes, (c) 30 minutes, (d) 12 hours (11.5 hours after injection stops). In each panel, the fault network geometry is shown in 2D view at the bottom, together with the fluid pressure (color scale on the right side). Slip velocity profiles are shown as curves above each fault, in logarithmic scale (vertical axis).</figcaption></figure>
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

Figure 7: Snapshots of slip velocity and fluid pressure as a function of distance from the injection point. The snapshot times are the same as in Fig. 6. In each panel, the extent of each fault is indicated by colored horizontal lines at the top, and their slip velocities by curves of the corresponding color. The black curve shows the maximum pressure across the fault network at each distance from the injection point.
<figure><img src="image_18.png" /></figure>

Figure 8: Distance from injection point of the slip velocity front defined by the threshold $V = 10^{-11}$ m/s (red curve). The transparent rectangles indicate the spatial extent of each fault that has reached the velocity threshold and its duration of slip at a velocity above the threshold. The red curves are selected fluid pressure contours; the darker red curve corresponds to the pressure front defined consistently with equation 17. The black dashed line indicates the end of the injection.

# 4 Conclusion

559

560

561

We presented a robust numerical framework, supported by an open-source software implementation, to simulate seismic and aseismic slip on non-planar single fault and planar
562

complex fault networks, specifically induced by fluid injection and diffusion along faults.

563

Our method seamlessly integrates a boundary element method, tailored to model slip on faults governed by rate-and-state friction, with a finite volume method for accurately capturing fluid diffusion along fault networks. To demonstrate the versatility and capabilities of our approach, we analyzed two representative scenarios: (1) fluid injection triggering slow slip on a primary, rough, rate-strengthening fault, which in turn initiates microseismicity on nearby secondary faults, and (2) fluid injection on one fault within a network of intersecting faults, leading to fluid diffusion and slip reactivation across the network.

564

565

566

567

568

569

570

These results highlight the intricate interplay between fluid diffusion and stress transfer by fluid-induced slow slip events. This complexity, evident in both scenarios, emphasizes the need to consider coupled hydro-mechanical processes when interpreting fault behavior and assessing seismic hazards.

571

572

573

574

This work provides a new simulator for exploring the interplay between fluid injection, fault mechanics, and seismic activity, which can offer valuable insights for understanding and mitigating induced seismicity in crustal reservoirs. While our model does not implement variable permeability yet, the current work paves the way to incorporate it in future studies, allowing for the full hydro-mechanically-coupled modeling in networks of geometrically complex faults. The permeability of faults is strongly influenced by normal stress (Barton et al., 1985; Esaki et al., 1999) and the main effect of complex fault geometry is to induce significant variations in normal stress during slip (Romanet et al., 2020, 2024). Future developments will incorporate variable fault permeability and poroelastic coupling to better capture the interplay between fault geometry and permeability changes, extending the model's applicability to more complex geological settings.

575

576

577

578

579

580

581

582

583

584
## Data Availability Statements

The software, input files and wrapper used to generate the figures, as well as the two videos of the main examples in this paper are available in a repository that can be accessed via the following DOI link: https://doi.org/10.5281/zenodo.15261403.

## Acknowledgments

P.R. thanks Barnaby Fryer for insights on the Finite Volume Method. P.R. would like to thank Harsha S. Bhat for his contribution in the early development of the code. This study was supported by European Research Council (ERC) Starting Grant 101040600 (HYQUAKE). J.P.A. was supported by the French government through the UCAJEDI Investments in the Future project (ANR-15-IDEX-01) managed by the National Research Agency (ANR). F.C. was supported by the Institut Universitaire de France. We kindly thank the two reviewers: Eric Dunham and an anonymous reviewer, as well as the Editor Shiqing Xu that helped us to improve the manuscript.
Almakari, M., Dublanchet, P., Chauris, H., and Pellet, F. (2019). Effect of the injection scenario on the rate and magnitude content of injection-induced seismicity: Case of a heterogeneous fault. *Journal of Geophysical Research: Solid Earth*, 124(8):8426–8448.

Ampuero, J.-P. and Rubin, A. M. (2008). Earthquake nucleation on rate and state faults–aging and slip laws. *J. Geophys. Res.*, 113(B1).

Atkinson, G. M., Eaton, D. W., and Igonin, N. (2020). Developments in understanding seismicity triggered by hydraulic fracturing. *Nature Reviews Earth & Environment*, 1(5):264–277.

Barton, N., Bandis, S., and Bakhtar, K. (1985). Strength, deformation and conductivity coupling of rock joints. In *International journal of rock mechanics and mining sciences & geomechanics abstracts*, volume 22, pages 121–140. Elsevier.

Berre, I., Doster, F., and Keilegavlen, E. (2019). Flow in fractured porous media: A review of conceptual models and discretization approaches. *Transport in Porous Media*, 130(1):215–236.

Bhattacharya, P. and Viesca, R. C. (2019). Fluid-induced aseismic fault slip outpaces pore-fluid migration. *Science*, 364(6439):464–468.

Bigaroni, N., Scuderi, M. M., Cappa, F., Guglielmi, Y., Nussbaum, C., Aldega, L., Pozzi, G., and Collettini, C. (2023). Frictional properties of opalinus clay: influence of humidity, normal stress and grain size on frictional stability. *Geophysical Journal International*, 233(1):211–228.

Börm, S., Grasedyck, L., and Hackbusch, W. (2003). Introduction to hierarchical matrices with applications. *Engineering analysis with boundary elements*, 27(5):405–422.
621

622

Bradley, A. M. (2014). Software for efficient static dislocation-traction calculations in fault simulators. *Seismol. Res. Lett.*, 85(6):1358-1365.

623

624

625

Candela, T., Brodsky, E. E., Marone, C., and Elsworth, D. (2014). Laboratory evidence for particle mobilization as a mechanism for permeability enhancement via dynamic stressing. *Earth and Planetary Science Letters*, 392:279-291.

626

627

Candela, T., Renard, F., Schmittbuhl, J., Bouchon, M., and Brodsky, E. E. (2011). Fault slip distribution and fault roughness. *Geophys. J. Int.*, 187(2):959-968.

628

629

630

Cappa, F., Guglielmi, Y., and De Barros, L. (2022a). Transient evolution of permeability and friction in a slowly slipping fault activated by fluid pressurization. *Nature communications*, 13(1):3039.

631

632

633

Cappa, F., Guglielmi, Y., Nussbaum, C., and Birkholzer, J. (2018). On the relationship between fault permeability increases, induced stress perturbation, and the growth of aseismic slip during fluid injection. *Geophysical Research Letters*, 45(20):11-012.

634

635

636

Cappa, F., Guglielmi, Y., Nussbaum, C., De Barros, L., and Birkholzer, J. (2022b). Fluid migration in low-permeability faults driven by decoupling of fault slip and opening. *Nature Geoscience*, 15(9):747-751.

637

638

639

Cappa, F., Scuderi, M. M., Collettini, C., Guglielmi, Y., and Avouac, J.-P. (2019). Stabilization of fault slip by fluid injection in the laboratory and in situ. *Science Advances*, 5(3):eaau4065.

640

641

Carslaw, H. S. and Jaeger, J. C. (1959). *Conduction of heat in solids*. Oxford Clarendon Press.

642

Chaillat, S., Desiderio, L., and Ciarlet, P. (2017). Theory and implementation of h-matrix
643

based iterative and direct solvers for helmholtz and elastodynamic oscillatory kernels. *J.*

644

*Comput. Phys.*, 351:165–186.

645

Ciardo, F. and Lecampion, B. (2023). Injection-induced aseismic slip in tight fractured

646

rocks. *Rock Mechanics and Rock Engineering*, 56(10):7027–7048.

647

Cochard, A. and Madariaga, R. (1996). Complexity of seismicity due to highly rate-

648

dependent friction. *J. Geophys. Res.*, 101(B11):25321–25336.

649

Cornet, F. H. (2012). The relationship between seismic and aseismic motions induced by

650

forced fluid injections. *Hydrogeology journal*, 20(8):1463.

651

Danré, P., Garagash, D., De Barros, L., Cappa, F., and Ampuero, J.-p. (2024). Control of

652

seismicity migration in earthquake swarms by injected fluid volume and aseismic crack

653

propagation. *Journal of Geophysical Research: Solid Earth*, 129(1):e2023JB027276.

654

De Barros, L., Danré, P., Garagash, D., Cappa, F., and Lengliné, O. (2024). Systematic ob-

655

servation of a seismic back-front during fluid injection in both natural and anthropogenic

656

earthquake swarms. *Earth and Planetary Science Letters*, 641:118849.

657

De Barros, L., Guglielmi, Y., Cappa, F., Nussbaum, C., and Birkholzer, J. (2023). Induced

658

microseismicity and tremor signatures illuminate different slip behaviours in a natural

659

shale fault reactivated by a fluid pressure stimulation (mont terri). *Geophysical Journal*

660

*International*, 235(1):531–541.

661

Dieterich, J. (1994). A constitutive law for rate of earthquake production and its application

662

to earthquake clustering. *J. Geophys. Res.*, 99(B2):2601–2618.

663

Dieterich, J. H., Richards-Dinger, K. B., and Kroll, K. A. (2015). Modeling injection-

664

induced seismicity with the physics-based earthquake simulator rsqsim. *Seismological*

665

*Research Letters*, 86(4):1102–1109.
666

Dublanchet, P. (2019). Fluid driven shear cracks on a strengthening rate-and-state frictional fault. *J. Mech. Phys. Solids*, 132:103672.

667

668

Duboeuf, L., De Barros, L., Cappa, F., Guglielmi, Y., Deschamps, A., and Seguy, S. (2017). Aseismic motions drive a sparse seismicity during fluid injections into a fractured zone in a carbonate reservoir. *Journal of Geophysical Research: Solid Earth*, 122(10):8285-8304.

669

670

671

Duboeuf, L., De Barros, L., Kakurina, M., Guglielmi, Y., Cappa, F., and Valley, B. (2021). Aseismic deformations perturb the stress state and trigger induced seismicity during injection experiments. *Geophysical Journal International*, 224(2):1464-1475.

672

673

674

Dunham, E. M. (2024). Fluid-driven aseismic fault slip with permeability enhancement and dilatancy. *Philosophical Transactions A*, 382(2276):20230255.

675

676

Dunham, E. M., Belanger, D., Cong, L., and Kozdon, J. E. (2011). Earthquake ruptures with strongly rate-weakening friction and off-fault plasticity, part 2: Nonplanar faults. *Bull. Seism. Soc. Am.*, 101(5):2308-2322.

677

678

679

Elkhoury, J. E., Brodsky, E. E., and Agnew, D. C. (2006). Seismic waves increase permeability. *Nature*, 441(7097):1135-1138.

680

681

Ellsworth, W. L. (2013). Injection-induced earthquakes. *Science*, 341.

682

Ellsworth, W. L., Giardini, D., Townend, J., Ge, S., and Shimamoto, T. (2019). Triggering of the pohang, korea, earthquake (m w 5.5) by enhanced geothermal system stimulation. *Seismological Research Letters*, 90(5):1844-1858.

683

684

685

Esaki, T., Du, S., Mitani, Y., Ikusada, K., and Jing, L. (1999). Development of a shear-flow test apparatus and determination of coupled properties for a single rock joint. *International journal of rock mechanics and mining sciences*, 36(5):641-650.

686

687
688

689

Evans, D. M. (1966). The denver area earthquakes and the rocky mountain arsenal disposal well. *The Mountain Geologist*.

690

691

692

Eyre, T. S., Eaton, D. W., Garagash, D. I., Zecevic, M., Venieri, M., Weir, R., and Lawton, D. C. (2019). The role of aseismic slip in hydraulic fracturing-induced seismicity. *Science advances*, 5(8):eaav7172.

693

694

695

Eyre, T. S., Samsonov, S., Feng, W., Kao, H., and Eaton, D. W. (2022). Insar data reveal that the largest hydraulic fracturing-induced earthquake in canada, to date, is a slow-slip event. *Scientific reports*, 12(1):2043.

696

697

Faoro, I., Elsworth, D., and Marone, C. (2012). Permeability evolution during dynamic stressing of dual permeability media. *J. Geophys. Res.*, 117:B01310.

698

699

Garagash, D. I. (2021). Fracture mechanics of rate-and-state faults and fluid injection induced slip. *Philosophical Transactions of the Royal Society A*, 379(2196):20200129.

700

701

Garagash, D. I. and Germanovich, L. N. (2012). Nucleation and arrest of dynamic slip on a pressurized fault. *J. Geophys. Res.*, 117(B10310).

702

703

704

Goebel, T., Weingarten, M., Chen, X., Haffener, J., and Brodsky, E. (2017). The 2016 mw5.1 fairview, okLAhoma earthquakes: Evidence for long-range poroelastic triggering at, 40 km from fluid disposal wells. *Earth Planet. Sc. Lett.*, 472:50-61.

705

706

707

708

Grigoli, F., Cesca, S., Priolo, E., Rinaldi, A. P., Clinton, J. F., Stabile, T. A., Dost, B., Fernandez, M. G., Wiemer, S., and Dahm, T. (2017). Current challenges in monitoring, discrimination, and management of induced seismicity related to underground industrial activities: A european perspective. *Reviews of Geophysics*, 55(2):310-340.

709

710

Grigoli, F., Cesca, S., Rinaldi, A. P., Manconi, A., Lopez-Comino, J. A., Clinton, J., Westaway, R., Cauzzi, C., Dahm, T., and Wiemer, S. (2018). The november 2017 m w
711

712

5.5 pohang earthquake: A possible case of induced seismicity in south korea. *Science*, 360(6392):1003-1006.

713

714

Guglielmi, Y., Cappa, F., Avouac, J.-P., Henry, P., and Elsworth, D. (2015). Seismicity triggered by fluid injection-induced aseismic slip. *Science*, 348(6240):1224-1226.

715

716

717

718

Guglielmi, Y., Nussbaum, C., Jeanne, P., Rutqvist, J., Cappa, F., and Birkholzer, J. (2020). Complexity of fault rupture and fluid leakage in shale: insights from a controlled fault activation experiment. *Journal of Geophysical Research: Solid Earth*, 125(2):e2019JB017781.

719

720

Hackbusch, W. (1999). A sparse matrix arithmetic based on matrices. part i: Introduction to matrices. *Computing*, 62(2):89-108.

721

722

723

Hamilton, D. H. and Meehan, R. L. (1971). Ground rupture in the baldwin hills: Injection of fluids into the ground for oil recovery and waste disposal triggers surface faulting. *Science*, 172(3981):333-344.

724

725

726

727

Heimisson, E. R., Liu, S., Lapusta, N., and Rudnicki, J. (2022). A spectral boundary-integral method for faults and fractures in a poroelastic solid: Simulations of a rate-and-state fault with dilatancy, compaction, and fluid injection. *Journal of Geophysical Research: Solid Earth*, 127(9):e2022JB024185.

728

729

730

731

Hicks, T., Pine, R., Willis-Richards, J., Xu, S., Jupe, A., and Rodrigues, N. (1996). A hydro-thermo-mechanical numerical model for hdr geothermal reservoir evaluation. In *International journal of rock mechanics and mining sciences & geomechanics abstracts*, volume 33, pages 499-511. Elsevier.

732

733

734

Im, K. and Avouac, J.-P. (2024). Quake-dfn: A software for simulating sequences of induced earthquakes in a discrete fault network. *Bulletin of the Seismological Society of America*, 114(5):2341-2358.
735
736

Im, K. and Avouac, J.-P. (2025). Maximum magnitude of induced earthquakes in rate and state friction framework. *Seismological Research Letters*, 96(3):1654-1664.

737
738

Im, K., Elsworth, D., and Fang, Y. (2018). The influence of preslip sealing on the permeability evolution of fractures and faults. *Geophys. Res. Lett.*, 45(1):166-175.

739
740

Im, K., Elsworth, D., and Wang, C. (2019). Cyclic permeability evolution during repose then reactivation of fractures and faults. *J. Geophys. Res.*, 124(5):4492-4506.

741
742

Karimi-Fard, M., Durlofsky, L. J., and Aziz, K. (2004). An efficient discrete-fracture model applicable for general-purpose reservoir simulators. *SPE journal*, 9(02):227-236.

743
744
745

Karvounis, D. and Wiemer, S. (2022). A discrete fracture hybrid model for forecasting diffusion-induced seismicity and power generation in enhanced geothermal systems. *Geophysical Journal International*, 230(1):84-113.

746
747
748

Katsman, R. and Ben-Avraham, Z. (2025). Modeling reservoir-induced seismicity: A dynamic poro-visco-elasto-plastic earthquake simulator with spontaneous dilatant coseismic rupture. *Earth and Space Science*, 12(2):e2024EA004029.

749
750
751

Keranen, K. M., Savage, H. M., Abers, G. A., and Cochran, E. S. (2013). Potentially induced earthquakes in okLAhoma, usa: Links between wastewater injection and the 2011 mw 5.7 earthquake sequence. *Geology*, 41(6):699-702.

752
753

Keranen, K. M. and Weingarten, M. (2018). Induced seismicity. *Ann. Rev. Earth Planet. Sci.*, (0).

754
755
756

Kim, K.-H., Ree, J.-H., Kim, Y., Kim, S., Kang, S. Y., and Seo, W. (2018). Assessing whether the 2017 m w 5.4 pohang earthquake in south korea was an induced event. *Science*, 360(6392):1007-1009.
757

Kohl, T. and Mégel, T. (2007). Predictive modeling of reservoir response to hydraulic stimulations at the european eggs site soultz-sous-forêts. *International Journal of Rock Mechanics and Mining Sciences*, 44(8):1118-1131.

758

759

760

Lai, G., Ge, H., Xue, L., Brodsky, E. E., Huang, F., and Wang, W. (2014). Tidal response variation and recovery following the wenchuan earthquake from water level data of multiple wells in the nearfield. *Tectonophysics*, 619:115-122.

761

762

763

Lambert, V. R., Erickson, B. A., Jiang, J., Dunham, E. M., Kim, T., Ampuero, J.-P., Ando, R., Cappa, F., Dublanchet, P., Elbanna, A., et al. (2025). Community-driven code comparisons for simulations of fluid-induced aseismic slip. *Journal of Geophysical Research: Solid Earth*, 130(4):e2024JB030601.

764

765

766

767

Larochelle, S., Lapusta, N., Ampuero, J.-P., and Cappa, F. (2021). Constraining fault friction and stability with fluid-injection field experiments. *Geophysical Research Letters*, 48(10):e2020GL091188.

768

769

770

Lee, H. and Cho, T. (2002). Hydraulic characteristics of rough fractures in linear flow under normal and shear load. *Rock mechanics and rock engineering*, 35:299-318.

771

772

Lee, K.-K., Ellsworth, W. L., Giardini, D., Townend, J., Ge, S., Shimamoto, T., Yeo, I.-W., Kang, T.-S., Rhie, J., Sheen, D.-H., et al. (2019). Managing injection-induced seismic risks. *Science*, 364(6442):730-732.

773

774

775

Lengliné, O., Boubacar, M., and Schmittbuhl, J. (2017). Seismicity related to the hydraulic stimulation of grt1, rittershoffen, france. *Geophysical Journal International*, 208(3):1704-1715.

776

777

778

Magnani, M. B., Blanpied, M. L., DeShon, H. R., and Hornbach, M. J. (2017). Discriminating between natural versus induced seismicity from long-term deformation history of intraplate faults. *Science advances*, 3(11):e1701593.

779

780
781

Manga, M., Beresnev, I., Brodsky, E. E., Elkhoury, J. E., Elsworth, D., Inge Britsen, S. E., Mays, D. C., and Wang, C. Y. (2012). Changes in permeability caused by transient stresses: Field observations, experiments, and mechanisms. *Rev. Geophys.*, 50(2):RG2004.

782

783

784

785

786

Marone, C. (1998). Laboratory-derived friction laws and their application to seismic faulting. *Ann. Rev. Earth Planet. Sci.*, 26(1):643-696.

787

788

789

Maurer, J., Dunham, E. M., and Segall, P. (2020). Role of fluid injection on earthquake size in dynamic rupture simulations on rough faults. *Geophysical Research Letters*, 47(13):e2020GL088377.

790

791

Maxwell, S., Zhang, F., and Damjanac, B. (2015). Geomechanical modeling of induced seismicity resulting from hydraulic fracturing. *The Leading Edge*, 34(6):678-683.

792

793

794

McClure, M. W. and Horne, R. N. (2011). Investigation of injection-induced seismicity using a coupled fluid flow and rate/state friction model. *Geophysics*, 76(6):WC181-WC198.

795

796

Noda, H. (2022). Dynamic earthquake sequence simulation with an sbiem accounting for interseismic poroelastic rebound. *Earth, Planets and Space*, 74(1):1-15.

797

798

799

800

Norbeck, J. and Horne, R. (2015). Injection-triggered seismicity: An investigation of porothermoelastic effects using a rate-and-state earthquake model. In *Proc., Fourtieth Workshop on Geothermal Reservoir Engineering, Stanford, California, USA*, pages 26-28.

801

802

803

Norbeck, J. and Horne, R. (2016). Evidence for a transient hydromechanical and frictional faulting response during the 2011 mw 5.6 prague, okLAhoma earthquake sequence. *Journal of Geophysical Research: Solid Earth*, 121(12):8688-8705.
804

Norbeck, J. H., McClure, M. W., Lo, J. W., and Horne, R. N. (2016). An embedded fracture modeling framework for simulation of hydraulic fracturing and shear stimulation. *Computational Geosciences*, 20(1):1–18.

805

806

807

Ohtani, M., Hirahara, K., Takahashi, Y., Hori, T., Hyodo, M., Nakashima, H., and Iwashita, T. (2011). Fast computation of quasi-dynamic earthquake cycle simulation with hierarchical matrices. *Procedia Comput. Sci.*, 4:1456–1465.

808

809

810

Orellana, L. F., Scuderi, M. M., Collettini, C., and Violay, M. (2018). Do scaly clays control seismicity on faulted shale rocks? *Earth and Planetary Science Letters*, 488:59–67.

811

812

Ozawa, S., Ida, A., Hoshino, T., and Ando, R. (2023). Large-scale earthquake sequence simulations on 3-d non-planar faults using the boundary element method accelerated by lattice h-matrices. *Geophysical Journal International*, 232(3):1471–1481.

813

814

815

Ozawa, S., Yang, Y., and Dunham, E. M. (2024). Fault-valve instability: A mechanism for slow slip events. *Journal of Geophysical Research: Solid Earth*, 129(10):e2024JB029165.

816

817

Passelègue, F. X., Brantut, N., and Mitchell, T. M. (2018). Fault reactivation by fluid injection: Controls from stress state and injection rate. *Geophysical Research Letters*, 45(23):12–837.

818

819

820

Patankar, S. (2018). *Numerical heat transfer and fluid flow*. CRC press.

821

Perfettini, H. and Ampuero, J.-P. (2008). Dynamics of a velocity strengthening fault region: Implications for slow earthquakes and postseismic slip. *J. Geophys. Res.*, 113(B9):B09411.

822

823

824

Perfettini, H. and Avouac, J.-P. (2007). Modeling afterslip and aftershocks following the 1992 landers earthquake. *Journal of Geophysical Research: Solid Earth*, 112(B7).

825
826

Power, W. L. and Tullis, T. E. (1991). Euclidean and fractal models for the description of rock surface roughness. *J. Geophys. Res.*, 96(B1):415–424.

827

828

Renard, F., Voisin, C., Marsan, D., and Schmittbuhl, J. (2006). High resolution 3d laser scanner measurements of a strike-slip fault quantify its morphological anisotropy at all scales. *Geophys. Res. Lett.*, 33(4).

829

830

831

Rice, J. R. (1993). Spatio-temporal complexity of slip on a fault. *J. Geophys. Res.*, 98(B6):9885–9907.

832

833

Romanet, P., Bhat, H. S., Jolivet, R., and Madariaga, R. (2018). Fast and slow slip events emerge due to fault geometrical complexity. *Geophys. Res. Lett.*, 45(10):4809–4819.

834

835

Romanet, P., Saito, T., and Fukuyama, E. (2024). The mechanics of static non-planar faults in infinitesimal strain theory. *Geophys. J. Int.*, page ggae337.

836

837

Romanet, P., Sato, D. S., and Ando, R. (2020). Curvature, a mechanical link between the geometrical complexities of a fault: application to bends, kinks and rough faults. *Geophys. J. Int.*

838

839

840

Rubin, A. and Ampuero, J.-P. (2005). Earthquake nucleation on (aging) rate and state faults. *J. Geophys. Res.*, 110:B11312.

841

842

Sáez, A. and Lecampion, B. (2024). Fluid-driven slow slip and earthquake nucleation on a slip-weakening circular fault. *Journal of the Mechanics and Physics of Solids*, 183:105506.

843

844

Sáez, A., Lecampion, B., Bhattacharya, P., and Viesca, R. C. (2022). Three-dimensional fluid-driven stable frictional ruptures. *Journal of the Mechanics and Physics of Solids*, 160:104754.

845

846

847

Sáez, A., Passelègue, F., and Lecampion, B. (2025). Maximum size and magnitude of injection-induced slow slip events. *Science Advances*, 11(19):eadq0662.

848
849

850

851

Sato, D. S. and Ando, R. (2021). A log-linear time algorithm for the elastodynamic boundary integral equation method. *Engineering Analysis with Boundary Elements*, 133:407-450.

852

853

854

Schoenball, M., Dorbath, L., Gaucher, E., Wellmann, J. F., and Kohl, T. (2014). Change of stress regime during geothermal reservoir stimulation. *Geophysical Research Letters*, 41(4):1163-1170.

855

856

857

Scotti, O. and Cornet, F. (1994). In situ evidence for fluid-induced aseismic slip events along fault zones. In *International journal of rock mechanics and mining sciences & geomechanics abstracts*, volume 31, pages 347-358. Elsevier.

858

859

860

Scuderi, M. M. and Collettini, C. (2016). The role of fluid pressure in induced vs. triggered seismicity: Insights from rock deformation experiments on carbonates. *Scientific reports*, 6(1):24852.

861

862

Shapiro, S. A., Huenges, E., and Borm, G. (1997). Estimating the crust permeability from fluid-injection-induced seismic emission at the ktb site. *Geophys. J. Int.*, 131(2).

863

864

865

Shapiro, S. A., Rothert, E., Rath, V., and Rindschwentner, J. (2002). Characterization of fluid transport properties of reservoirs using induced microseismicity. *Geophysics*, 67(1):212-220.

866

867

868

Sumy, D. F., Cochran, E. S., Keranen, K. M., Wei, M., and Abers, G. A. (2014). Observations of static coulomb stress triggering of the november 2011 m5. 7 okLAhoma earthquake sequence. *Journal of Geophysical Research: Solid Earth*, 119(3):1904-1923.

869

870

Tada, T. and Yamashita, T. (1997). Non-hypersingular boundary integral equations for two dimensional non-planar crack analysis. *Geophys. J. Int.*, 130(2):269-282.
871

Tan, X. and Lui, S. K. Y. (2025). Fluid-induced earthquake nucleation on aging rate-and-state faults: Influence of hydraulic diffusivity and injection rate under different nucleation regimes.

872

873

874

Townend, J. and Zoback, M. (2000). How faulting keeps the crust strong. *Geology*, 28(5):399.

875

876

Uphoff, C., May, D. A., and Gabriel, A.-A. (2023). A discontinuous galerkin method for sequences of earthquakes and aseismic slip on multiple faults using unstructured curvilinear grids. *Geophysical Journal International*, 233(1):586–626.

877

878

879

Viesca, R. C. (2021). Self-similar fault slip in response to fluid injection. *J. Fluid Mech.*, 928:A29.

880

881

Wang, T. A. and Dunham, E. M. (2022). Hindcasting injection-induced aseismic slip and microseismicity at the cooper basin enhanced geothermal systems project. *Scientific Reports*, 12(1):19481.

882

883

884

Wei, S., Avouac, J.-P., Hudnut, K. W., Donnellan, A., Parker, J. W., Graves, R. W., Helmberger, D., Fielding, E., Liu, Z., Cappa, F., et al. (2015). The 2012 brawley swarm triggered by injection-induced aseismic slip. *Earth and Planetary Science Letters*, 422:115–125.

885

886

887

888

Woo, J.-U., Kim, M., Sheen, D.-H., Kang, T.-S., Rhie, J., Grigoli, F., Ellsworth, W. L., and Giardini, D. (2019). An in-depth seismological analysis revealing a causal link between the 2017 mw 5.5 pohang earthquake and eggs project. *Journal of Geophysical Research: Solid Earth*, 124(12):13060–13078.

889

890

891

892

Xue, L., Li, H.-B., Brodsky, E. E., Xu, Z.-Q., Kano, Y., Wang, H., Mori, J. J., Si, J.-L., Pei, J.-L., Zhang, W., et al. (2013). Continuous permeability measurements record healing inside the wenchuan earthquake fault zone. *Science*, 340(6140):1555–1559.

893

894
895

896

897

Yang, Y. and Dunham, E. M. (2023). Influence of creep compaction and dilatancy on earthquake sequences and slow slip. *Journal of Geophysical Research: Solid Earth*, 128(4):e2022JB025969.

898

899

Yerkes, R. F. and Castle, R. O. (1976). Seismicity and faulting attributable to fluid extraction. *Engineering Geology*, 10(2-4):151-167.

900

901

902

Zhu, W., Allison, K. L., Dunham, E. M., and Yang, Y. (2020). Fault valving and pore pressure evolution in simulations of earthquake sequences and aseismic slip. *Nature communications*, 11(1):4833.
903 A Stress transfer kernels for BEM

The kernels relating slip to elastic shear stress $\tau_{el}(\mathbf{x})$ and normal stress $\sigma_{el}(\mathbf{x})$ are (Tada and Yamashita, 1997; Romanet et al., 2018, 2020):

$$ K_{\text{grad}}^{t}(\mathbf{x}, \mathbf{y}) = [4n_{1}(\mathbf{x})n_{2}(\mathbf{x})\gamma_{1}(\mathbf{x}, \mathbf{y})\gamma_{2}(\mathbf{x}, \mathbf{y}) + (n_{2}^{2}(\mathbf{x}) - n_{1}^{2}(\mathbf{x}))\gamma_{2}^{2}(\mathbf{x}, \mathbf{y}) - \gamma_{1}^{2}(\mathbf{x}, \mathbf{y})] \\ \times \left(n_{2}(\mathbf{y})\frac{\gamma_{1}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})} - n_{1}(\mathbf{y})\frac{\gamma_{2}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})}\right), $$

$$ K_{\text{curv}}^{t}(\mathbf{x}, \mathbf{y}) = [4n_{1}(\mathbf{x})n_{2}(\mathbf{x})\gamma_{1}(\mathbf{x}, \mathbf{y})\gamma_{2}(\mathbf{x}, \mathbf{y}) + (n_{2}^{2}(\mathbf{x}) - n_{1}^{2}(\mathbf{x}))\gamma_{2}^{2}(\mathbf{x}, \mathbf{y}) - \gamma_{1}^{2}(\mathbf{x}, \mathbf{y})] \\ \times \left(-n_{1}(\mathbf{y})\frac{\gamma_{1}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})} - n_{2}(\mathbf{y})\frac{\gamma_{2}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})}\right), $$

$$ K_{\text{grad}}^{n}(\mathbf{x}, \mathbf{y}) = - \left( n_{1}(\mathbf{y})\frac{\gamma_{1}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})} + n_{2}(\mathbf{y})\frac{\gamma_{2}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})} \right. \\ \left. + [2n_{1}(\mathbf{x})n_{2}(\mathbf{x})(\gamma_{2}^{2}(\mathbf{x}, \mathbf{y}) - \gamma_{1}^{2}(\mathbf{x}, \mathbf{y})) - 2\gamma_{1}(\mathbf{x}, \mathbf{y})\gamma_{2}(\mathbf{x}, \mathbf{y})(n_{2}^{2}(\mathbf{x}) - n_{1}^{2}(\mathbf{x}))] \right. \\ \left. \times \left(n_{2}(\mathbf{y})\frac{\gamma_{1}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})} - n_{1}(\mathbf{y})\frac{\gamma_{2}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})}\right) \right), $$

$$ K_{\text{curv}}^{n}(\mathbf{x}, \mathbf{y}) = - \left( n_{2}(\mathbf{y})\frac{\gamma_{1}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})} - n_{1}(\mathbf{y})\frac{\gamma_{2}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})} \right. \\ + [2n_{1}(\mathbf{x})n_{2}(\mathbf{x})(\gamma_{2}^{2}(\mathbf{x}, \mathbf{y}) - \gamma_{1}^{2}(\mathbf{x}, \mathbf{y})) - 2\gamma_{1}(\mathbf{x}, \mathbf{y})\gamma_{2}(\mathbf{x}, \mathbf{y})(n_{2}^{2}(\mathbf{x}) - n_{1}^{2}(\mathbf{x}))] \\ \left. \times \left(-n_{1}(\mathbf{y})\frac{\gamma_{1}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})} - n_{2}(\mathbf{y})\frac{\gamma_{2}(\mathbf{x}, \mathbf{y})}{r(\mathbf{x}, \mathbf{y})}\right) \right) $$

where $\gamma_{1}(\mathbf{x}, \mathbf{y}) = \frac{x_{1}-y_{1}}{r(\mathbf{x}, \mathbf{y})}$, $\gamma_{2}(\mathbf{x}, \mathbf{y}) = \frac{x_{2}-y_{2}}{r(\mathbf{x}, \mathbf{y})}$, and $r(\mathbf{x}, \mathbf{y}) = \sqrt{(x_{1}-y_{1})^{2} + (x_{2}-y_{2})^{2}}$ is the distance between the two points $\mathbf{x}$ and $\mathbf{y}$. The convention that compressive stresses are positive has been adopted here, contrary to classic boundary integral methods.

909 B SCEC SEAS benchmark problem BP06

We benchmarked the code using the SCEC SEAS benchmark problem BP6 of fluid injection at constant flow rate on a planar velocity-strengthening rate-and-state fault with the aging law (Lambert et al., 2025). The parameter values are listed in Tab. 3 and a more complete
913 description of the benchmark problem can be found on the SCEC SEAS website: https://strike.scec.org/cvws/seas/benchmark_descriptions.html. Contrary to the two 914 examples described in the main text, the benchmark problem involves antiplane slip (mode 915 III). On planar faults, the governing equations of inplane and antiplane modes are closely 916 related: the mode II equations are the same as the mode III ones after replacing the shear 917 modulus $\mu$ by $\mu/(1 - \nu)$. Thus, the mode III benchmark considered here serves also as a 918 benchmark for mode II on planar faults. 919
<table><caption>Table 3: Table of parameters for the benchmark SCEC SEAS BP6 with the aging law (Lambert et al., 2025).</caption><tbody><tr><td colspan="3">Frictional fault parameters</td></tr><tr><td>a</td><td>0.007</td><td>Direct effect</td></tr><tr><td>b</td><td>0.005</td><td>Evolution effect</td></tr><tr><td>D<sub>c</sub></td><td>4 mm</td><td>Critical Distance</td></tr><tr><td>f<sub>0</sub></td><td>0.6</td><td>Reference friction</td></tr><tr><td>V<sub>0</sub></td><td>10<sup>-6</sup> m s<sup>-1</sup></td><td>Reference velocity</td></tr><tr><td colspan="3">Mechanical properties for rocks</td></tr><tr><td>μ</td><td>32.04 GPa</td><td>Shear modulus</td></tr><tr><td>c<sub>s</sub></td><td>3.464 km s<sup>-1</sup></td><td>Shear wave speed</td></tr><tr><td colspan="3">Fault hydraulic and fluid properties</td></tr><tr><td>α</td><td>0.1 m<sup>2</sup> s<sup>-1</sup></td><td>Hydraulic diffusivity</td></tr><tr><td>k</td><td>10<sup>-13</sup> m<sup>2</sup></td><td>Permeability</td></tr><tr><td>β<sub>f</sub> + β<sub>φ</sub></td><td>10<sup>-8</sup> Pa<sup>-1</sup></td><td>Pore and fluid compressibility</td></tr><tr><td>φ</td><td>0.1</td><td>Porosity</td></tr><tr><td>η</td><td>10<sup>-3</sup> Pa.s</td><td>Dynamic viscosity</td></tr><tr><td>q</td><td>1.25 × 10<sup>-6</sup> m s<sup>-1</sup></td><td>Flow rate of injection</td></tr><tr><td>t<sub>inj</sub></td><td>0 s</td><td>Beginning of injection</td></tr><tr><td>Δt<sub>inj</sub></td><td>100 days</td><td>Duration of injection</td></tr><tr><td colspan="3">Initial values for mechanical conditions</td></tr><tr><td>V<sub>ini</sub></td><td>10<sup>-12</sup> m s<sup>-1</sup></td><td>Initial slip velocity</td></tr><tr><td>θ<sub>ini</sub></td><td>4.0956 × 10<sup>10</sup> s</td><td>Initial state</td></tr><tr><td>σ<sub>load</sub></td><td>50 MPa</td><td>Normal stress</td></tr><tr><td colspan="3">Lengthscales</td></tr><tr><td>ds</td><td>10 m</td><td>Grid size</td></tr><tr><td>L<sub>b</sub></td><td>513 m</td><td>Process zone size</td></tr></tbody></table>
<figure><img src="image_19.png" /><figcaption>Figure 9: Snapshots of the SCEC SEAS BP6 benchmark simulation. Selected times, relative to the beginning of injection, are: (a) 50 days, (b) 100 days, (c) 150 days (this is when injection stops)), and (d) 2 years.</figcaption></figure>

Fig. 9 shows different snapshots of the fluid diffusion along the fault, and the slow slip it triggered. In Fig. 10 we compare the results of our code to those of two other codes: the Discontinuous Galerkin Finite Element Method (DGFEM) software TANDEM (Uphoff et al. (2023), with grid size of 100 m) and the spectral Boundary Element Method (SBEM) (Dublanchet (2019), with grid size of 9.76 m). The results of the three codes closely agree, both for the slip velocity (Figs. 10.a and 10.c) and for the state variable (Figs. 10.b and 10.d).

We verified our FVM implementation by comparing its results against the analytical solution for pore pressure diffusion in the injection scenario of the BP06 benchmark (see supporting information of Lambert et al. (2025)). The results of our FVM closely match the analytical solution (Fig. 11). Their relative difference (normalized $L_2$ norm shown in Fig. 12.a for different grid sizes) is largest ($\sim 1\%$) at the beginning of injection because the pressure gradient is very strong at the injection point, and somewhat grid-dependent. It then decreases over time. Long after the end of injection, the pressure profile is flat and
the relative error becomes very small (<< 0.1 %). At long times, the relative error converges approximately linearly (figure 12.b) as a function of decreasing grid size, as expected for the first-order backward Euler scheme.

<figure><img src="image_20.png" /></figure>

Figure 10: (a) Slip velocity and (b) state evolution at the center of the fault ($z = 0$ km). (c) Slip velocity and (d) state evolution at $z = 2.5$ km. The results of our BEM-FVM code are compared to those of a discontinuous Galerkin method (DGM) (Uphoff et al., 2023) and a spectral boundary element method (SBEM) (Dublanchet, 2019). The vertical dashed lines indicate the end of injection.
<figure><img src="image_21.png" /><figcaption>Figure 11: Comparison of pressure profiles computed by the Finite Volume Method (FVM) and the analytical solution of the fluid diffusion equation under injection at constant flow rate at $x = 0$.</figcaption></figure>
<figure><img src="image_22.png" /></figure>

Figure 12: a) Error of the FVM relative to the analytical solution for different fault discretizations, as a function of time. As diffusion progresses, pressure becomes smoother and the relative error decreases. b) Error of the FVM relative to the analytical solution as a function of gridsize. The error seems to decrease linearly with time except for the roughest discretisation.

## C Benchmark for crossing faults

We conducted benchmark simulations on two crossing faults with varying grid size to verify the convergence of the numerical solution. The two faults are 10 m long and cross at their center with an angle $\phi$ (Fig. 13). Injection at constant flow rate is applied during a finite time interval at a point on fault 1, and the two fault have a no flow condition at their edges. Model parameters are given in Tab. 4.
<figure><img src="image_23.png" /><figcaption>Figure 13: Fault geometry for the benchmark of two crossing faults.</figcaption></figure>
<table><caption>Table 4: Parameter values for the benchmark on crossing faults.</caption><thead><tr><th colspan="3">Fault frictional properties</th></tr></thead><tbody><tr><td>a</td><td>0.005</td><td>Direct effect</td></tr><tr><td>b</td><td>0.004</td><td>Evolution effect</td></tr><tr><td>D<sub>c</sub></td><td>1 µm</td><td>Critical slip distance</td></tr><tr><td>f<sub>0</sub></td><td>0.6</td><td>Reference friction</td></tr><tr><td>V<sub>0</sub></td><td>1 µm s<sup>-1</sup></td><td>Reference velocity</td></tr><tr><th colspan="3">Elastic properties of rock</th></tr><tr><td>μ</td><td>32.040 GPa</td><td>Shear modulus</td></tr><tr><td>ν</td><td>0.1541</td><td>Poisson's ratio</td></tr><tr><th colspan="3">Fault hydraulic and fluid properties</th></tr><tr><td>k</td><td>1 × 10<sup>-15</sup> m<sup>2</sup></td><td>Permeability</td></tr><tr><td>β<sub>f</sub></td><td>9 × 10<sup>-9</sup> Pa<sup>-1</sup></td><td>Compressibility of the fluid</td></tr><tr><td>β<sub>φ</sub></td><td>1 × 10<sup>-9</sup> Pa<sup>-1</sup></td><td>Compressibility of the pores</td></tr><tr><td>φ</td><td>0.1</td><td>Porosity</td></tr><tr><td>η</td><td>1 × 10<sup>-3</sup> Pa.s</td><td>Dynamic viscosity</td></tr><tr><td>q</td><td>2 × 10<sup>-6</sup> m s<sup>-1</sup></td><td>Injection flow rate</td></tr><tr><td>t<sub>inj</sub></td><td>3 years</td><td>Beginning of injection</td></tr><tr><td>Δt<sub>inj</sub></td><td>2000 s</td><td>Duration of injection</td></tr><tr><th colspan="3">Initial values for mechanical variables</th></tr><tr><td>V<sub>ini</sub></td><td>1 × 10<sup>-12</sup> m s<sup>-1</sup></td><td>Initial slip velocity</td></tr><tr><td>θ<sub>ini</sub></td><td>10<sup>6</sup> s</td><td>Initial state</td></tr><tr><td>σ<sub>load</sub></td><td>5 MPa</td><td>Normal stress</td></tr><tr><th colspan="3">Length scales</th></tr><tr><td>L<sub>b</sub></td><td>1.67 m</td><td>Process zone size</td></tr></tbody></table>
Given the no flow-condition for the edges of the network, the final pressure of the simulations can be calculated following:

$$P_{\text{final}} = \frac{Q * \Delta t_{\text{inj}}}{\phi(\beta_f + \beta_\phi)L_{\text{faults}}}, \qquad (19)$$

where $\Delta t_{\text{inj}}$ is the total duration of injection, and $L_{\text{faults}}$ is the cumulative length of the faults in the network. For this two-crossing fault benchmark, the final pressure is therefore $P_{\text{final}} = 0.2$ MPa.

Figures 14 to 16 show the spatial distributions of slip velocity and fluid pressure at two selected times after injection. Each figure is for a different angle between the faults: $\phi = 30^{\circ}$, $60^{\circ}$ and $90^{\circ}$, respectively. Simulation results with different grid sizes are shown: $ds = 0.005, 0.01$ and $0.05$ m. The coarsest grid size is ~ 33 times larger than the theoretical process zone size $L_b$. Because the time step is adaptive, the results shown for each simulation are at the output time closest to the time indicated above each panel. The results from all grid sizes overlap. This indicates that the numerical solution converges with decreasing grid size. The final pressures from the simulations (Figs. 14.f, 15.f and 16.f) is also converging correctly towards the theoretical final pressure $P_{\text{final}} = 0.2$ MPa.
<figure><img src="image_24.png" /><figcaption>Figure 14: Results of the benchmark problem with two faults crossing at an angle $\phi = 30^\circ$, for different grid sizes (see legend at the bottom): slip velocity (a, c, and d) and pore-pressure (b,d, and f) along fault 1 (blue) and fault 2 (red), at three different times indicated above each plot.</figcaption></figure>
a) 0003 years, 000 days, 00 hours, 10 minutes, 00.000 seconds

<figure><img src="image_25.png" /></figure>

b) 0003 years, 000 days, 00 hours, 10 minutes, 00.000 seconds

<figure><img src="image_26.png" /></figure>

c) 0003 years, 000 days, 00 hours, 33 minutes, 20.000 seconds

<figure><img src="image_27.png" /></figure>

d) 0003 years, 000 days, 00 hours, 33 minutes, 20.000 seconds

<figure><img src="image_28.png" /></figure>

e) 0003 years, 001 days, 00 hours, 00 minutes, 00.000 seconds

<figure><img src="image_29.png" /></figure>

f) 0003 years, 001 days, 00 hours, 00 minutes, 00.000 seconds

<figure><img src="image_30.png" /><figcaption>Figure 15: Same as Fig. 14 but for $\phi = 60^{\circ}$.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure 16: Same as Fig. 14 but for $\phi = 90^{\circ}$.</figcaption></figure>

D Additional simulation results on triggered microseismicity

To assess the effect of the secondary faults on the main fault in the example shown in section 3.1, we conducted the same simulation but without the secondary faults. The
results of the two simulations are compared in Fig. 17. When a secondary fault ruptures,

it has an effect (due to static stress transfer) on the slip velocity of the main fault (Figs.

17.a and 17.c). However, this effect vanishes very quickly: the slip velocities of the main

fault in the two simulations are almost identical soon after the rupture of a secondary fault

(Figs. 17.b and 17.d).
<figure><img src="image_32.png" /><figcaption>Figure 17: Comparison of slip velocities on the main fault obtained in two simulations with (blue) and without (red) secondary faults, at four different times (indicated above each plot). (a) One hour after injection, the effect of a secondary fault rupturing can be seen on the slip velocity of the main fault. (b) Five minutes later, this effect has vanished. (c and d) Another example, in which the effect of a secondary fault rupturing on the main fault has vanished after 14 minutes.</figcaption></figure>

To check the robustness of the seismicity features highlighted in section 3.1, we ran
967

additional simulations using the same parameter values, but changing the random seed for

968

the generation of the rough main fault and of the orientation of the secondary faults. The

969

resulting spatio-temporal distributions of seismicity are shown in Fig. 18. Fig. 18.a is

970

the same simulation as presented in section 3.1. Similar seismicity patterns are observed:

971

seismicity migrates away from the main fault during injection, and continues long after the

972

injection stopped.
<figure><img src="image_33.png" /><figcaption>Figure 18: Microseismicity migration, same as Fig. 4 but for 6 different random realizations of secondary fault orientations and main fault roughness. Each panel shows results from a single simulation.</figcaption></figure>