We use a 3D Distinct Element Code (3DEC, Itasca Consulting Group, Inc, 2016) to model a fluid injection in a permeable, slip-weakening fault separating two three-dimensional (3D) elastic, impervious blocks (Figure 1). We consider the distinct element method (Cundall, 1988) to simulate the hydromechanical interactions and slip along the fault. The blocks are meshed in tetrahedral zones linked by nodes. The fault is discretized in 0.7 m-sized cells. Using an explicit time-marching procedure, the distinct element method considers force-displacement relations updated at each time step to describe the interaction between the blocks. Newton’s laws (i.e., the differential equations of motion) and constitutive relations are used to calculate the translation and rotation of the blocks at each time step, and thus their position, velocity, displacement and applied nodal forces. The discontinuities forces are updated through constitutive laws before next time step.

The model is hydromechanically fully coupled, with hydraulic and frictional parameters that vary during injection. This affects the normal and shear displacement on the fault as well as the seismic and aseismic deformation during rupture. Seismicity is generated using an inherently discrete rupture model (Wynants-Morel et al., 2020). Given the values adopted in our model, the nucleation size of earthquakes for a slip-weakening fault (Uenishi and Rice, 2003) varies between 0.01 and 0.32 m, depending on the effective normal stress at rupture and the frictional fault properties (Tables 1, 2). It means that the nucleation phase of the earthquakes is smaller than the cell size, which prevents to compute it accurately, but the model allows to generate seismic events with rupture size of a m² up to few tens of m². We use an adaptive time step to be able to simulate both fast ruptures (i.e., seismic) and slow deformations (i.e., aseismic). The time step can decrease down to 10^–6 second during co-seismic ruptures and increase to 1 s during aseismic periods. We then estimate earthquake source properties using classical seismological assumptions detailed in Wynants-Morel et al. (2020).

During fluid injection into a fault, the pressure and fluid flow are related to the hydraulic aperture through the cubic law (Witherspoon et al., 1980). Before the injection, an initial hydraulic aperture (b_ho) is assigned to the fault within the initial background stress and pressure conditions. Variations of hydraulic aperture (Δb_h) from initial aperture come from two main processes: (1) variations of effective normal stress and (2) dilations while the fracture slips. Thus, they can be expressed as:

where Δσ_n’ is the increment in effective normal stress (total normal stress minus fluid pressure) (Pa), k_n is the normal stiffness (Pa/m), Δu_s (m) is the shear slip increment, and ψ is the dilation angle (°).

The permeability (k) is related to the square of the hydraulic aperture (Jaeger and Cook, 1984). Modeling studies show that the permeability on a fault affects the fault slip behavior (Cappa and Rutqvist, 2011; Yeo et al., 2020) but also the spatio-temporal distribution of seismic events (Shapiro et al., 1997). In addition, as the permeability increases with increasing pressure and fault slip, the evolution of permeability during injection can affect aseismic deformation developing outside the pressurized zone (Cappa et al., 2018) as well as induced seismicity at short and long distances from injection (Rinaldi and Nespoli, 2017). Thus, the study of the permeability and its variations are fundamental to understand the relationship between injection parameters and the seismic and aseismic slip released during fluid injection.

The existence of a strong coupling between permeability and shear failure is frequently observed in natural reservoirs and in-situ experiments (e.g., Guglielmi et al.,2015a,b; Zhang and Li, 2016). Previous studies showed that the dilation angle is adequate to model this coupling, that occurs during deep underground CO₂ injection (Cappa and Rutqvist, 2011), the stimulation of a geothermal reservoir (Fomin et al., 2004) or in-situ injections (Guglielmi et al., 2015b; Tsopela et al., 2019). This parameter is commonly used to model the simulation of the geomechanical and hydraulic response of natural fracture systems (Latham et al., 2013; Lei et al., 2014, 2016). It is thus a critical parameter to understand hydromechanical coupling during fault deformation.

The Mohr-Coulomb failure criterion (Jaeger and Cook, 1984) is given by:

where τ_f is the shear strength of the fault (Pa), σ_n is the normal stress (Pa), and P is the fluid pressure (Pa) acting on the fault. μ_S is the friction coefficient before failure. The fault is reactivated when the Mohr-Coulomb criterion is reached (i.e., when the shear stress and the shear strength of the faults are equal). Thus, slip begins and the slip velocity may increase. We assume here a null cohesive strength of the fault (c = 0).

The Shear Capacity Utilization (SCU), proposed in Buijze et al. (2019a), is a parameter that allows us to quantify the initial closeness to failure of the fault. It is defined as the ratio between initial shear stress (τ₀) and initial strength on the fault (τ_f0), before the injection starts:

with σ_{n_0} and P₀ the initial normal stress and pressure acting on the fault. A maximal SCU (100%) corresponds to a fault at failure that can slip without any perturbation.

The Mohr-Coulomb criterion can be reached by increasing either the fluid pressure (and, thus decreasing the effective stress and the shear strength) or the shear stress. We use the Shear Stress Contribution (SSC) to quantify the respective role of the fluid pressure and the shear stress in the rupture:

where τ0⁢(x¯),σn0⁢(x¯), and P0⁢(x¯) are, respectively, the shear stress, the normal stress and the fluid pressure measured at position x¯ before injection. △⁢τ⁢(x¯) is the shear stress variation between the initial state and the state for which rupture occurs at the location x¯. The SSC varies between 0 and 100%, respectively between a failure only triggered by fluid pressure changes and a failure only induced by shear stress perturbations.

During the slip periods, we consider that the friction coefficient only evolves with the shear slip, assuming a conventional linear slip-weakening model (Ida, 1972). In this law, a friction drop occurs as the friction coefficient decreases linearly over a critical slip distance (D_c) from a peak static value (μ_S) to a residual dynamic value (μ_D) (Figure 2A and Table 1). Both the friction drop (μ_S−μ_D) and D_c control the slip behavior (e.g., Ampuero et al., 2002; Mikumo et al., 2003; Uenishi and Rice, 2003; Dunham, 2007), and thus the seismic or aseismic nature of the fault response.

In our modeling, the use of a linear slip-weakening friction law reflects of a simplified approach of fault slip and earthquake rupture modeling. However, for simplicity, the slip-weakening friction is here preferred to the more sophisticated rate-and-state friction law as our model is used to show a series of few seismic events in a large aseismic rupture (Wynants-Morel et al., 2020).

Our objective is to study the effect of fault physical properties on the distribution of seismic events and the seismic-to-total partitioning during fluid injection in a single fault, at a typical depth of storage reservoirs (3 km). This fault, with a dip angle of 70°, is 106 m long and 100 m large, along dip and strike, respectively. Around the fault, the rock medium is elastic and impervious (Figure 1). In the parametric analysis, we consider a reference case from which the other simulations differ by varying one parameter (Table 2).

We apply constant principal stresses to the model (σ_zz = 81 MPa, σ_xx = σ_yy = 70 MPa), which vary with depth through a lithostatic gradient due to gravity (g = 9.81 m/s²). We then test different initial stress state, varying the SCU of the fault between 53 and 71%. In order to do so, we add a deviatoric stress component σ_xy that varies between 13.5 and 18.5 MPa. It allows us to keep the normal stress constant among the tests and change the shear stress only. Before injection, we apply an initial fluid pressure (P₀) of 30 MPa, which varies with depth with a hydrostatic gradient.

To activate the fault, we consider a fluid injection at the center of the fault (Figure 1). The applied flowrate linearly increases during the first hundred seconds of injection and it is then kept constant at a value of 0.09 m³/s (Figure 2B). This short-duration, high-rate injection leads to a fast pressurization of the fault at the injection point. Nonetheless, considering a single small injection point may not be pertinent in certain operational injection scenarios such as saltwater disposal where multiple wells are generally used. The simulations are stopped when the rupture over the fault reached a fixed distance from injection point (corresponding to 38 m and 76% of the length of the fault). Therefore, the simulation duration differs among tests, and lasts 190 s in our reference case (Figure 2B). However, considering an imposed rupture size implies that the behavior of the seismicity is well-represented at short distances from the injection, and therefore for very short duration lengths.

We set elastic properties for rocks to typical constant values for sedimentary material (Table 1). These values are standard in crustal reservoirs, as for instance in Oklahoma, a region of pronounced injection-induced seismicity (Barbour et al., 2017). We test different values of dilation angle (Table 2), which are consistent with previous simulations (Segall et al., 2010; McClure and Horne, 2011; Guglielmi et al., 2015b). Several values of fault initial permeability are also tested. Nonetheless, the considered values are very high, to explore the behavior of highly conductive channels that can be observed in subsurface reservoirs (Jeanne et al., 2012).

We consider a classical static friction coefficient for faults at crustal conditions (μ_S = 0.6; Byerlee, 1978). We then test different values of the dynamic friction coefficient μ_D in the range 0.375–0.475. The critical slip distance D_c also varies between 5 and 80 μm. Such values are consistent with laboratory observations (Marone, 1998; Rubino et al., 2017). No fault healing is considered as test durations are small in this study, and also because healing is not well-known during fluid injection.

In a simulation, the evolution of fluid pressure, fault opening, fault slip, stress, friction, and hydraulic aperture are monitored. Seismic ruptures are detected using a typical threshold on the slip velocity from Cochard and Madariaga (1994):

where c_S is the shear wave speed (m/s) and G is the rock shear modulus (Pa). Usual velocity threshold values from literature go from 0.1 mm/s to 0.1 m/s (McClure and Horne, 2011; Gischig, 2015; McClure, 2015). We adopt here a 1 mm/s threshold.

Neighboring grid cells that show subsequent seismic slip velocities during an overlapping period of time are then grouped to form seismic events when at least seven cells are regrouped (Wynants-Morel et al., 2020). The seismic event stops when the slip velocity of the last cell drops under the threshold. Non-seismic ruptures are considered aseismic. At the end of a simulation, the location, timing and moment magnitude of seismic events are estimated using the criteria described in Wynants-Morel et al. (2020).

First, we examine the common set of hydromechanical and seismological features between the 21 injection tests considered in our study. The first rupture occurred at the injection point, then the rupture starts expanding radially on the fault (Wynants-Morel et al., 2020). The resulting fault slip patch has an elliptic shape elongated in the maximal shear stress direction (Figure 3). Oppositely, the fluid diffuses with a circular shape in all cases (Wynants-Morel et al., 2020). The induced deformation is a combination of aseismic slip and seismic events. The maximum aseismic slip is located at the injection point (4.1–18.0 mm, depending on the simulation), and decreases progressively when getting closer to the rupture edge. Induced by the fluid injections, some seismic events occurred. In the simulations, a minimum threshold for the moment magnitude of seismic events has been fixed to –1.6 to have rupture on patches greater than 6 grid cells. This limit is equivalent to a detection threshold for induced seismicity recorded at reservoir-scale (e.g., Kwiatek et al., 2019). Thus, the distribution of seismic slip on the fault is sparse (Figure 4). We observe 7–73 events, depending on the test, with magnitude ranging from –1.6 to 0.3. Cumulated seismic slip is smaller than a few tens of microns (0.114–0.383 mm, depending on the simulation). It is 20–128 times smaller than the maximum aseismic slip (0.413–1,800 mm). Therefore, the calculated deformation is mostly aseismic, with at most 4% of the total moment that is seismic. This result is consistent with independent seismological observations in both in-situ decametric-scale and small-scale laboratory experiments (e.g., Goodfellow et al., 2015; De Barros et al., 2018). As the rupture is stopped at a fixed distance from injection, the ruptured area and final aseismic moment are very similar for all simulations.

To better understand the relationship between seismicity, fluid pressure and shear stress, we represent the seismicity in a distance-time diagram (or r-t plot) (Figure 5). The pressure front, defined by an overpressure of 5% compared to the hydrostatic level, follows a diffusion curve of the form r = √4πDt (Shapiro et al., 1997) where D is a diffusivity coefficient (m²/s). As the fault friction reduces with slip (i.e., slip-weakening), when the rupture occurs, the slipping patch has a residual shear stress lower than the initial one, and is surrounded by an area with increased stress. We therefore define the shear stress front as the peak of stress at the edge of the slipping area. Since the rupture zone has an elliptical shape, the distance from the shear stress front to the injection varies with the angles on the fault. Therefore, in a r-t plot, the shear stress front varies within a band that gets wider as the rupture grows with time (Figure 5). The rupture front is defined as the largest distance of the shear stress front. We observe that the seismicity always occurs on the shear stress front and does not follow the pressure front (Figure 5). The shear stress front may or may not outpace the pressure front. When it does, the seismic front accelerates, and tends to have a constant or increasing migration velocity. This migration velocity corresponds to the propagation velocity of the aseismic rupture. A pronounced acceleration of the seismic moment rate is observed at the same time. This feature is observed for large SCU and friction drop tests (Figures 5I,L). On the contrary, when the shear stress front stays behind the pressure front, the seismic migration velocity slows down, and shows a diffusive shape.

In this study, injection is flowrate-controlled. Thus, time and injected volume are interdependent (Figure 2B). Since the test duration evolves with the initial permeability, the friction drop and the SCU (Figures 5D–L), the injected volume needed to reach a same rupture size depends on these three fault properties. Moreover, since the rupture front at the end of the simulation is fixed, an increase in the test duration implies a decrease in the migration velocity of the shear stress and seismicity fronts.

The shear stress contribution required to induce failure (SSC, Eq. 4) is calculated to quantify the respective contribution of shear stress and fluid pressure in the triggering of seismic events. It increases with time and distance to injection (Figure 5). In 10 out of 21 simulations, the first 1–to–3 seismic events have a SSC under 50% during each test, which means they are mostly triggered by fluid overpressure. In the other simulations, the SSC is always between 50 and 99%, which corresponds to seismic events dominantly triggered by shear stress variations.

Variations in the SCU and friction drop lead to very similar behaviors (Figures 5G–L). First, increasing SCU or friction drops, leads to smaller fluid injected volume, and smaller test durations, without any changes in the pressure front shape (Table 3). For small SCU or friction drop, the shear stress front and thus the seismicity stay behind the pressure front. The seismicity migrates with small velocities (less than 0.1 m/s) which slow down with time (Figure 6F). For larger values, the shear stress front is first behind the pressure front, before accelerating and outpacing it. Therefore, the seismic migration velocity increases with increasing SCU or friction drop. The shear stress contribution (SSC, Eq. 4), the number of seismic events and the cumulative seismic moment at the end of injection also increase with increasing SCU and friction drop (Figures 6D–E, 7). Therefore, the seismic-to-total moment ratio also increases with increasing SCU (from 0.040 to 3.2%) and friction drop (from 0.034 to 4.2%).

Variations among tests with different dilation angle (Figures 5A–C, 6C–E, 7 and Table 2) are not observed in the spatio-temporal distribution of the seismicity, the shear stress front, the seismic moment, the seismic-to-total moment ratio, nor in the number of events. The main differences when varying the dilation angles concern the amplitude of the fluid pressure at the injection, and the shape of the pressure front. With a null dilation angle (i.e., no shear-induced dilation), the diffusive-like pressure front reaches 29 m at the end of injection (Figure 5A), while with a 8° dilation angle, it stops following a diffusion curve near 100 s and accelerates to reach 39 m (Figure 5C). Therefore, while the stress front outpaces the pressure front for low dilation angles, it stays behind for high dilation angles.

While an increase of initial fault permeability allows a faster diffusion of the fluid, it also reduces the migration velocity of the rupture (Figure 6F) and the seismicity (Table 3). Therefore, an increase of the initial permeability leads to an increase of the injected fluid volume and test duration to reach the same rupture size (Figures 5D–F). Moreover, permeability increases in the pressurized zone during the simulation. As shown in Figure 8A, the permeability variations increase at injection point with decreasing initial permeability (from 53 to 2,500% of its initial value) and increasing dilation angle (from 108 to 1,700%). Nonetheless, the seismicity and shear stress fronts still accelerate and outpace the pressure front during the injection simulation. Increasing initial permeability also induces a decrease in the final cumulative seismic moment released during the simulation (from 47.0 × 10⁸ N.m to 5.8 × 10⁸ N.m), as well as in the seismic-to-total moment ratio (from 1.8 to 0.38%). The initial permeability mainly acts on the fluid pressure at the injection, with a change of behavior that can be observed for k = 1.78 × 10^–9 m². For smaller values, the increase of pressure at the injection leads to a null effective normal stress before the end of the simulation (Figure 6B). Therefore, a maximum number of events is reached for this value (Figure 6E).

A factor of 10 increase in the critical slip distance D_c only affects the number of seismic events and the cumulative seismic moment which are, respectively, divided by 10 and 16 (Figures 6E, 7). Therefore, the seismic-to-total moment ratio is also strongly reduced when increasing critical slip distance. Neither the spatiotemporal repartition of the seismicity, the shear stress front nor the fluid pressure front are influenced by this parameter (Table 3).

Here we examine how the operational parameters that can be monitored (i.e., injected fluid volume, seismicity) during reservoir stimulations depend on the fault properties.

First, for a same rupture size, the injected fluid volume is strongly modified by the friction drop, the SCU and the initial permeability (Figure 6), highlighting the importance of the hydromechanical coupling on the fault ruptured area. It is divided by 21.4 for a SCU multiplied by 1.3 and by 10.2 for a friction drop multiplied by 1.8. It is also multiplied by 12.75 for a 10-fold increase of the initial permeability. No significant injection volume variations are observed with dilation angle and critical slip distance simulations.

The pressure front at the end of the injections can be modeled using a single hydraulic diffusivity (D) equals to 0.37 m²/s for the simulations with varying friction drop, SCU and critical slip distance (Figure 6A). Therefore, these three fault properties do not influence the hydrological behavior of the fault. A saturation of the pressure front is nonetheless observed for the smallest SCU test (SCU = 53%), as the pressure front reaches the edges of the model. On the contrary, the diffusivity measured on the pressure front obviously varies with the initial permeability, from 0.28 to 0.68 m²/s. A non-null dilation angle also changes the shape of the pressure front, which cannot be modeled by a linear diffusion law anymore (see section “Influence of Fault Physical Properties on the Fault Sismo-Hydromechanical Responses”). Besides, both initial permeability and dilation angles strongly modify fluid pressure at the injection (Figure 6B). Indeed, after a sharp increase of the injection pressure leading to the first rupture, the pressure at the injection slowly increases when no dilation is considered, but it decreases with dilation angles larger than 2° (Figure 8B). For smaller dilation angles, the pressure at injection point decreases later, after the rupture of a larger zone around injection. Moreover, the maximum overpressure reached during the simulations decreases with an increasing dilation angle (38.4–17.2 MPa for a dilation angle rising from 0° to 8°) and with an increasing initial permeability (46.3–25.9 MPa from k = 5.62 × 10^–10 m² to k = 5.62 × 10^–9 m²). The SCU and friction drop also influence the fluid pressure at the injection, but to a lesser extent.

A population of seismic events is induced during each simulation. The number of events is mostly impacted by the critical slip distance, the SCU and the friction drop (Figure 6E). It is divided by 10.4 for an 8-fold increase of the critical slip distance. This number is also multiplied by 5 and 3.7 for a SCU and friction drop, respectively, multiplied by 1.3 and 1.8. The initial permeability also acts on the number of events, but due to the pressure threshold reached in high permeability tests, no particular evolution can be shown regarding this fault property (see section “Influence of Fault Physical Properties on the Fault Sismo-Hydromechanical Responses”). No noticeable variations of the number of events are observed with dilation angle simulations.

The cumulative seismic moment is mostly influenced by the friction drop and the SCU, and, to a lesser extent, by the initial permeability and the critical slip distance (Figure 7). It is multiplied by 41.6 and 76.2 for a SCU and friction drop multiplied, respectively, by 1.3 and 1.8. It is also divided by 8.3 and 15.6 for an 8-fold increase of the critical slip distance and a 10-fold increase of the initial permeability, respectively. No significant seismic moment variations are observed with dilation angle simulations. As the total deformation is similar for all simulations, the seismic-to-total moment partitioning has the same sensitivity as the seismic moment, and is mainly modified by the friction drop and the SCU. The seismic-to-total moment ratio is indeed divided by 80 and 120 within the explored range of these two parameters. The critical slip distance and the initial permeability have a smaller influence (factor 16 and factor 5, respectively).

The shear stress contribution (SSC, Eq. 4) quantifies how much the shear stress (and conversely, the fluid pressure) contributes to the seismic rupture. Friction drop and SCU strongly influence the SSC. SSC decreases quickly for most events when friction drop and SCU decreases, meaning that the fluid pressure is playing a larger role in inducing seismic slips over the course of injection (Figure 6D). Besides, we also note an influence of the dilation angle on the SSC, which varies from 98 to 83% for dilation angles between 0° and 8°.

Therefore, the spatiotemporal distribution of induced seismicity as well as the partitioning between seismic and aseismic moment released during injection are deeply influenced by the initial proximity of stress to failure of the fault and by its friction drop with slip, and, to a lesser extent, by the initial permeability of the fault and the critical slip distance D_c. These behaviors are consistent with independent field observations (De Barros et al., 2016, 2018; Lund Snee and Zoback, 2016; Hearn et al., 2018).

During a fluid injection, the released seismic moment is usually related to the injected fluid volume (McGarr, 2014; Galis et al., 2017; McGarr and Barbour, 2018). The 21 simulations we have performed can be seen as the response of 21 different reservoirs. We showed that the volume depends on three fault parameters (i.e., the SCU, the friction drop, the initial permeability; Figure 7). The seismic moment also depends on these three parameters, as well as on the critical slip distance. In a moment-volume diagram (Figure 7), the dependence we observed between these two parameters is counter-intuitive, and not in agreement with classical laws. The seismic moment decreases with the injected volume when three parameters (the initial permeability, the SCU and the friction drop) are modified. Therefore, the effect of the considered parameters on the fault hydraulic, mechanical and seismic responses has to be considered in order to anticipate the released seismic moment from the injected fluid volume.

We find in our models that an increase in the friction drop or the SCU leads to a faster migration of the rupture and of the seismicity (Figures 5J–L) and to an increase of both the magnitude and the number of seismic events (Figure 7). These fault properties act on the background stress drop, defined as the difference between the initial stress state and the residual shear stress μ_d × σ_n after an aseismic or seismic rupture (e.g., Galis et al., 2017). The background stress drop represents the shear stress exceeding the dynamic strength of the fault before injection. It corresponds to an excess or a lack of stress, which can either lead or slow down the rupture, depending of the sign of the background stress drop. If the friction drop or the SCU are high, the background stress drop is positive. Therefore, the aseismic rupture, as well as the shear stress front, accelerates. The increase of the seismicity migration velocity with the SCU is observed in previous modeling studies considering a pressurized fault driven by a slip-weakening (Bhattacharya and Viesca, 2019; Wynants-Morel et al., 2020) or a rate-and-state friction (Gischig, 2015; Dublanchet, 2019a). Moreover, the acceleration of the seismicity migration in high SCU cases is also mentioned in previous studies (Dublanchet, 2019a; Wynants-Morel et al., 2020).

Both the number and magnitudes of induced earthquakes are larger with higher SCU and friction drop. This increase of the seismic moment with SCU and friction drop was observed in previous works (Gischig, 2015; Rutqvist et al., 2015; Dublanchet, 2019a; Wynants-Morel et al., 2020). Moreover, the SCU acts on the triggering of seismic events and on their spatial distribution (Keranen et al., 2014; Lund Snee and Zoback, 2016; Jin and Zoback,2018a,b; Buijze et al., 2019b). Indeed, for instance in Oklahoma (US), earthquakes are triggered even for stress and fluid pressure perturbations less than 3 MPa, because of the near-critical stress state of the faults (Lund Snee and Zoback, 2016; Schoenball and Ellsworth, 2017).

The initial permeability and the dilation angle both modify the hydraulic response of the fault. As shown by Eq. (1), while the evolution of permeability is related to the initial opening of the fault and the change in effective stress, the dilation angle only acts on the opening of the fault during the slip episodes. The initial permeability mainly acts on the migration velocity of seismic events, but oppositely to what is expected from poroelastic diffusion (Shapiro et al., 1997). An increase of initial permeability indeed decreases the seismic migration velocity (Figure 6F). It can be explained by the lower level of fluid pressure reached on the fault at a given time for higher initial permeability. As the injection is flowrate-controlled in this study (Figure 2B), the fluid pressure increases more slowly in a higher permeability simulation. Moreover, the SSC, and thus the pressure contribution, is not impacted by the initial permeability (Figure 6D). Consequently, the duration needed to reach a fixed level of fluid pressure at the maximum rupture distance increases with the permeability, which means a slower rupture migration velocity. This is in contradiction with Chang and Segall (2016) and Norbeck and Horne (2016), who showed that the rupture is faster with a higher permeability. Nonetheless, shear stress variations from aseismic slip is not considered in their studies, while it contributes to the majority of the rupture propagation in our model (Figure 6D; Wynants-Morel et al., 2020). It also supports the importance of full hydromechanical coupling, which is not taken into account in Chang and Segall (2016). Our study also shows that the cumulative seismic moment, and therefore the seismic-to-total ratio, increases with decreasing fault permeability (Figure 7). Such behavior has been observed in reservoir injection sites, for instance in Oklahoma (Shah and Keller, 2017; Hearn et al., 2018).

We observe that the Shear Stress Contribution to reach rupture (SSC) decreases with an increase of the dilation angle. This is due to a higher opening of the fault before rupture due to the slipping of close regions; therefore, fluid pressure can diffuse more easily and reaches higher values just before rupture. Thus, the pressure contribution to reach rupture is higher, which means a smaller SSC. Moreover, as the dilation angle does not influence the rupture itself, dilation angles have no effect on the seismic moment and the spatio-temporal evolution of the seismic front (Figure 7). This does not agree to what has been found in previous models (Segall et al., 2010; McClure and Horne, 2011; Ciardo and Lecampion, 2019), in which a damping of the seismicity is observed. However, these models do not consider shear stress variations from previous rupture as a mechanism to reach rupture. This highlights the importance of shear stress variations in the seismic triggering mechanism.

The cumulated seismic moment decreases with increasing critical slip distance D_c (Figure 7). Moreover, a high value of this frictional parameter means that the friction drops on a larger slip distance (Figure 2A). Therefore, the friction weakening rate (i.e., the ratio between the friction drop and the critical slip distance) is small, which impedes slip acceleration and reduces the possibilities of unstable, seismic slip (Hillers et al., 2009). Thus, a sparse seismicity is observed, leading to a small released seismic moment. A similar result is shown in Cueto-Felgueroso et al. (2017) for an injection in a rate-and-state fault. Rutqvist et al. (2015) observed a decrease of the maximal seismic moment for an increasing D_c, but also an increase of the number of events. Nonetheless, they do not consider aseismic rupture in their model; most of the events observed by Rutqvist et al. (2015) for high critical slip distance may be considered aseismic in our model.

We now investigate the contribution and consequences of our modeling results on the seismic moment vs. injected fluid volume relationship (Figure 7) commonly used to estimate the maximum earthquake magnitude induced by a fluid injection (Galis et al., 2017; McGarr and Barbour, 2018).

McGarr (1976) first considered the seismic moment (M0seis) released during a fluid injection to be related at first order to the injected fluid volume (V) through the shear modulus (G) of the surrounding rock. Then, this relationship has been updated in McGarr and Barbour (2018) to incorporate the aseismic moment (M0aseis) occurring during the injection:

This relationship stands for an upper bound that cannot be exceeded during a fluid injection.

Our results show that four fault parameters (i.e., the initial permeability, the SCU, the friction drop and the critical slip distance) modify the moment-volume relationship. Indeed, the seismic moment decreases with the injected volume when three parameters (the initial permeability, the SCU and the friction drop) are modified. The sensitivity to the critical slip distance of the fault leads to similar injected volume but different seismic moment (Figure 7). Assuming that these 21 tests are responses from different and independent reservoirs, the mean trend that shows up from this test is that the released moment decreases with injected volume. This is in total disagreement with the McGarr’s relations. Because our simulations are stopped for similar total moment, the latter does not follow either the eq. 6. However, even if the total moment was corrected to fit such a law, the seismic moment would not show an increasing dependency with the volume. Indeed, we here observed that seismic and total moment do not follow the same trend, i.e., that the seismic-to-total ratio globally decreases here with increasing injected fluid volume. It means that the fault physical properties have a significant impact on the moment-volume relationships by modifying the seismic-to-total ratio, at least within the small range of volume and seismic moment we have here explored. To better anticipate the released seismicity, the aseismic component of the deformation should be considered (De Barros et al., 2019). This motivates a reformulation of the available moment-volume models since most of them do not include the aseismic deformation that occurs during fluid injection. Nonetheless, we want to precise that the seismic moment does follow a power law with the injected volume in each individual test, as observed in the literature (Buijze et al., 2015; Bentz et al., 2020).

The relationship given by McGarr and Barbour (2018) (Eq. 6) is considering the total moment, that is the seismic plus aseismic deformation. However, it is exceeded for small values of initial permeability and it will probably be exceeded for values of SCU and friction drop larger than those adopted in this study (Figure 7). This is in agreement with reservoir stimulation cases in which the observed moment is far above the expected value obtained from Eq. (6) (Atkinson et al., 2016; Ellsworth et al., 2019). For instance, Chang et al. (2020) consider that the low permeability of the basement near the geothermal field of Pohang, South Korea, favored the triggering of the 2017 post-injection M_w 5.4 event. According to Atkinson et al. (2016), the diffusion of fluids from previous near injections brought a large zone in a critical state of stress close to the location of hydraulic fracturing exploitations in the WCSB. This may explain that 6 M_w > 3.5 events induced by these exploitations do not follow the prevision of the McGarr’s law (Atkinson et al., 2016). Therefore, seismic moment that largely exceeds prediction from the fluid volume may exist in case of low permeability or highly critical stress state.

As neither the seismic nor the total moment released during injection can be directly anticipated from the injected fluid volume only, we look for another measurable parameter that can be considered as a proxy of the released seismic moment. We observe that no detectable seismic event is observed close to injection point (Figure 5), which is consistent with De Barros et al. (2016) and Duboeuf et al. (2017) observations in in-situ controlled decametric scale experiments. We approximate the size of the aseismic slip patch around injection point as the distance between the injection point and the first seismic event with a magnitude above a magnitude threshold M_t = –1.15. Following the approach of De Barros et al. (2019), we observe that the induced moment of events with magnitude M_w ≥ M_t decreases with the distance to first seismic event for the 16 simulations (Figure 9) with seismic events with magnitude above –1.15. A scaling relating the seismic moment to the power –3 of the size of the aseismic patch around injection give a poor fit with the data at first order. Therefore, this parameter seems to be a possible probe to estimate the amount of aseismic deformation. Nonetheless, the power coefficient is lower to the coefficient proposed by De Barros et al. (2019), equal to –1.5. This difference can be explained by pressure or stress effects, as De Barros et al. (2019) find a different coefficient (equal to –1) by including the critical fluid pressure into the moment/distance relation. Therefore, such parameter, in combination with other parameters, might be useful to correct the moment-volume relationship for the fault properties and the aseismic deformation. Nonetheless, this parameter stays dependent to the knowing of the injection location, which may not be available in many deep injection cases.