Impact of injection rate ramp-up on nucleation and arrest of dynamic fault slip

Fluid injection into underground formations reactivates preexisting geological discontinuities such as faults or fractures. In this work, we investigate the impact of injection rate ramp-up present in many standard injection protocols on the nucleation and potential arrest of dynamic slip along a planar pressurized fault. We assume a linear increasing function of injection rate with time, up to a given time tc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_c$$\end{document} after which a maximum value Qm\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q_m$$\end{document} is achieved. Under the assumption of negligible shear-induced dilatancy and impermeable host medium, we solve numerically the coupled hydro-mechanical model and explore the different slip regimes identified via scaling analysis. We show that in the limit when fluid diffusion time scale tw\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_w$$\end{document} is much larger than the ramp-up time scale tc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_c$$\end{document}, slip on an ultimately stable fault is essentially driven by pressurization at constant rate. Vice versa, in the limit when tc/tw≫1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_c/t_w \gg 1$$\end{document}, the pressurization rate, quantified by the dimensionless ratio QmtwtcQ∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dfrac{Q_m t_w}{t_c Q^*}$$\end{document} with Q∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q^*$$\end{document} being a characteristic injection rate scale, does impact both nucleation time and arrest distance of dynamic slip. Indeed, for a given initial fault loading condition and frictional weakening property, lower pressurization rates delay the nucleation of a finite-sized dynamic event and increase the corresponding run-out distance approximately proportional to ∝QmtwtcQ∗-0.472\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto \left( \dfrac{Q_m t_w}{t_c Q^*}\right) ^{-0.472}$$\end{document}. On critically stressed faults, instead, the ramp-up of injection rate activates quasi-static slip which quickly turn into a run-away dynamic rupture. Its nucleation time decreases non-linearly with increasing value of QmtwtcQ∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dfrac{Q_m t_w}{t_c Q^*}$$\end{document} and it may precede (or not) the one associated with fault pressurization at constant rate only.


Introduction
Anthropogenic fluid injection into underground formations is a common operation in many industrial applications. In the context of deep geothermal energy extraction, for instance, fluid is injected into targeted deep fault/fracture zones in order to enhance reservoir permeability and hence fluid circulation between injection and production wells (Giardini 2009;Deichmann and Giardini 2009). Among other applications that involve injection of fluid into subsurface there are hydraulic fracturing for oil and gas extraction from hydrocarbon reservoirs and wastewater disposal using deep wells (Warpinski and Teufel 1987;Horton 2012 Although these engineering techniques are widely used, the injection of fluids into underground formations alters the local equilibrium of the Earth's crust, inducing micro-seismicity and, in some cases, large earthquakes (Horton 2012;Ellsworth 2013;Kim 2013;Keranen et al. 2014;Weingarten et al. 2015). Significant earthquakes have been directly correlated with the injection activities, among these the Pohang earthquake of 2017 in South Korea with M w ¼ 5:5 (Kim et al. 2018;Grigoli et al. 2018;Yeo et al. 2020), the four M w ¼ 3 events in Basel, Switzerland, between 2006and 2007(Deichmann and Giardini 2009Goertz-Allmann et al. 2011), an event of M w ¼ 3:5 in the city of St. Gallen, Switzerland, back in 2013 (Diehl et al. 2014;Edwards et al. 2015;Diehl et al. 2017) and the M w ¼ 5:7 earthquake near Prague, Oklahoma, in 2011 (Keranen et al. 2013;Sumy et al. 2014).
Many numerical and theoretical models have been developed in order to investigate the impact of operational design parameters, such as injection pressure or injection rate, on fault slip activation and earthquakes nucleation upon fluid injection. Most of them are based on a Rate-and State-dependent friction model and, therefore, are well suited to explain features of earthquake cycles and seismicity rates. For example, Dempsey and Riffault (2019) used a pressure diffusion model coupled to R&S friction model to show that a reduction in injection rate may lead to a decrease in the seismicity rate in Oklahoma (USA). A similar result has been obtained by Lagenbruch and Zoback (2016) using instead a statistical model calibrated over many injection induced-earthquakes in Oklahoma. Using a poroelastic model incorporating R&S friction, Barbour et al. (2017) observed that in Oklahoma a variable injection rate may lead to a larger seismicity rate increases compared to the one under constant injection rate (for an equivalent injected volume). Chang et al. (2018), instead, studied the effect of injection rate variation on seismicity rate post shut-in and they showed that a gradual reduction of injection rate minimises the postinjection seismicity rate. Using a Dietrich-Ruina heterogeneous 2D fault, Almakari et al. (2019) investigated the effect of injection scenario not only in terms of seismicity rate, but also in terms of magnitude content. They showed that the total seismic moment increases with both maximum pressure and pressure rate and that the total number of induced seismic events is controlled by the maximum pressure. A recent study of Rudnicki and Zhan (2020) on a springblock model shows that larger pressurization rates stabilize fault slip events due to rate and state friction.
The role of injection design parameters on fault slip behaviour has been extensively investigated also in many laboratory experiments. Among others, Wang et al. (2020) showed that fault slip propagation is manly governed by fluid pressurization rate rather than injection pressure. French et al. (2016), instead, observed that fluid pressurization is less effective than mechanical changes in the fault normal stress at initiating accelerated slip events. The effect of fluid pressure oscillations on fault slip stability has been investigated by Noël et al. (2019) via a triaxial laboratory experiment. They showed that perturbations caused by pore fluid oscillations promote seismic slip and that seismic activity along a fault increases for increasing oscillation's amplitudes.
Despite the numerous studies on the effects of injection parameters on earthquakes nucleation and occurrence along faults, the impact of pressurization rate on the onset of dynamic fault slip remains still elusive. Garagash and Germanovich (2012) investigated extensively the conditions of nucleation and arrest of dynamic fault slip on a frictional weakening pressurized 2D fault. This work was then extended by Galis et al. (2017) to planar 3D faults subjected to point source injection and axisymmetric diffusion. Their generic findings, however, are valid only for two types of fault pressurizations, constant over-pressure and constant injection rate, which are over-simplification of many injection protocols commonly used in industrial applications. In numerous fault reactivation experiments, for instance related to deep geothermal energy exploitation, the injection protocol consists of one (or more) stimulation cycle, where the controlled injection rate or injection pressure increases in time (typically with stair-like increments), up to reach a steady state regime, followed then by a shut-in phase (Hofmann et al. 2018). In the hydraulic stimulation experiments conducted in Grimsel Test Site, Switzerland, back in 2017 for example, the injection protocol consisted of 4 injection cycles in which either injection pressure or injection flow rate was increased in a stepwise manner, before reaching a plateau and shut-in phase (Villiger et al. 2020). A similar trend of injection rate was also used during the hydraulic stimulation of Enhanced Geothermal System in the city of Basel, Switzerland (2006) (Häring et al. 2008) and in the more recent EGS project in Pohang, South Korea (Hofmann et al. 2018), to cite a few examples.
In this contribution, we extend the model of Garagash and Germanovich (2012) to account for an initial ramp-up of injection rate in time and investigate its effect on dynamic fault slip nucleation and arrest. We approximate the step-wise increase of injection rate adopted in standard injection protocols using a simple linear increasing function with time, followed by a maximum plateau (the shut-in phase is out of context and thus is not considered in this work). This choice represents a good approximation when the time scale of each increment is much larger that the corresponding one of each step (and thus a linear increasing function is a reasonable approximation). We solve the two-dimensional hydro-mechanical problem numerically and verify the results with theoretical predictions. The goal of this study is to evaluate the effects of the pressurization rate, quantified by the injection rate ramp-up variation, on the nucleation and arrest of a seismic rupture on a frictional weakening planar fault.

Fault model
We consider a planar fault embedded in an isotropic, homogeneous and unbounded elastic medium under plane-strain conditions (see Fig. 1). The fault is subjected to an ambient pore-pressure p o and a uniform far-field stress state that resolved on the fault plane result in an effective normal r 0 o ¼ r n À p o and shear s o stress component. Such a uniform ambient stress state, typical of a limited fault extent compared to the background in-situ gradient, is perturbed via a point source injection of volumetric flow rate Q(t) L 2 =T ½ directly in the middle of the fault (specifically at x ¼ 0). In order to investigate the effect of injection rate ramp-up on fault slip stability, we consider a linear increase of injection flow rate in time, followed by a plateau after a given time t c (the shut-in phase is out of scope here, hence it will not be considered-see Fig. 1). This parametrisation represents an approximation of many standard fault injection protocols used in hydro-shearing stimulation of fractured reservoirs, in which the design (constant) value of injection rate is reached upon stair-like increments.
Prior fluid injection, we assume that the fault is in static equilibrium with the uniform in-situ stress state (locked status). This tacitly assumes that there is no Fig. 1 Plane-strain fault model subjected to a farfield stress state and fluid injection. The conductive planar fault is embedded into an homogeneous, isotropic and linear-elastic medium characterized by a negligible hydraulic diffusivity. The thin fault gouge unit, therefore, accommodates accelerating slip due to friction weakening condition and pore-pressure diffusion effect of remote plates loading that would cause steady movements with slow energy release (creep) (Chen and Bürgmann 2017). The ambient equilibrium is violated throughout pressurization, during which pore-pressure perturbation diffuses along the fault plane and activates a symmetric shear crack of length 2a. This scenario of mechanically weak, fault unit accommodating slip and pore fluid flow may be representative of a highly permeable thick unit in a damage zone of a fault zone, whose hydraulic conductivity is much larger than the one of the host rock around it (and thus impermeable host medium is a reasonable assumption).
Before presenting the governing equations, it must be noted that we assume normal stresses positive in compression and shear stresses positive for clockwise rotation. Furthermore, thermal effects, dilatancy/compaction during shear deformations and poroelastic stress changes in the surrounding medium are neglected.
2.1 Governing equations for quasi-static slip development driven by pore-pressure diffusion 2.1.1 Static equilibrium and constitutive law for frictional slip We consider the activation and propagation of a symmetric shear crack of length 2a driven by pore fluid flow inside the conductive fault plane. The shear stress s at a given time t and position x on the slip surface is linearly related to slip d (or shear displacement discontinuity) within the slipping region via the following quasi-static elastic equilibrium equation where a þ and a À are respectively the positive and negative positions of the crack tips, s o is the uniform background shear stress, Young's modulus (with G and m being shear modulus and Poisson's ratio, respectively), odðf; tÞ=of is the shear dislocation density along the crack that must satisfy the following condition lim f!a AE ðtÞ odðf; tÞ of ffiffiffiffiffiffiffiffiffiffi ffi a Ç f p ¼ 0 ð2Þ in order to remove stress singularity at crack tips (Uenishi and Rice 2003), and 1 f À x is the non-local elastic kernel. For merely planar frictional problems that do not account for dilatancy or compaction during crack propagation, the elastic kernel affects only the shear stress distribution along the fault plane, while the uniform total normal stress r n remains constant. Inside the sliding region, the force balance requires that the shear stresses must equal the available frictional resistance which we assume here to obey the Mohr-Coulomb yielding criterion (without cohesion), accounting for a slip weakening of friction coefficient where ðr n À pðx; tÞÞ ¼ r 0 n ðx; tÞ is the local effective normal stress, strictly function of pore-pressure evolution p(x, t) inside the fault, and f ðdÞ is the slip weakening friction coefficient The frictional resistance, therefore, weakens linearly with shear deformations d, from a peak value associated with peak friction coefficient f p , to a residual value, after sufficiently large slip d r (such that the friction coefficient drops to its residual value f r ). Although more complicated frictional laws can be considered in (3), we use here the simple piece-wise linear weakening of friction coefficient with slip, which can be shown to be a good approximation of the phenomenological rate-and state-friction law at large slip rates and for small values of a/b (Uenishi and Rice 2003; Rubin and Ampuero 2005), or when Df p =b ) 1, with Df p being the peak change of friction coefficient from the steady-state value (Garagash 2021). Since the total weakening Df p on rate-and state-faults depends on the slip front speed, the condition Df p =b ) 1 holds when the front speed varies appreciably (see Fig. 4a in (Garagash 2021)) and thus the linear slip weakening friction law can realistically approximate the rate-and state-friction law.

Pore-pressure diffusion
Under the assumption of negligible fluid leak-off in the surrounding elastic medium, fluid flow is confined within the fault gouge unit characterised by a constant hydraulic aperture w h . Upon injection of volumetric flow rate Q(t) at x ¼ 0, the diffusion of pore fluid overpressure pðx; tÞ ¼ pðx; tÞ À p o is governed by the width-averaged fluid mass conservation equation where c f M À1 TL 2 ½ is a parameter that combines pore fluid compressibility and pore space expansivity and v is the gap-averaged fluid flow velocity given by the where k f L 2 ½ is the longitudinal fault permeability and l ¼ 12l 0 MT À1 L À1 ½ is a viscosity parameter (with l 0 as dynamic viscosity of the fluid).
The point injection boundary condition requires that where the volumetric flow rate Q(t) L 2 T À1 ½ increases linearly with time t, up to reach a plateau after a given time t c , i.e.
Since we neglect dilatation/compaction of fault gouge unit during shear crack propagation, fault permeability k f and hydraulic aperture w h remain constant throughout pressurization. Equations (5) and (6), therefore, reduce to the well-known parabolic diffusion equation that govern over-pressure diffusion inside the fault conduit characterized by a constant hydraulic diffu- Using specific boundary and initial conditions representative of the particular injection scenario and its time history, equation (9) can be solved analytically for the spatial and temporal evolution of pore-fluid over-pressure along the fault plane pðx; tÞ (see ''Appendix'' A for full details). During the ramp-up phase of injection rate, i.e. for t t c , the over-pressure evolution in function of time t and the normalized coordinate n ¼ x ffiffiffiffiffiffi ffi 4at p is given by where Erfc is the complementary error function. Notice that the analytical solution (10) is expressed as a product two independent functions, DPðtÞ and WðnÞ, which identify respectively the maximum overpressure evolution at injection point and its instantaneous spatial distribution. The pore-fluid evolution after the ramp-up phase, instead, is governed by the following equation (Cole et al. 2011 where ; t c denotes the over-pressure distribution at time t c and Gðx À x 0 ; t À t 0 Þ is the fundamental heat conduction solution valid for an infinite one-dimensional body subjected to an instantaneous point source (also called Green's function) (Carslaw and Jaeger 1959). We revert the reader to ''Appendix'' A for its analytical expression. It is worth mentioning that the solution (11) is valid for t [ t c and it takes into account the whole injection history during the rampup phase. If t c vanishes, then we recover the analytical solution of pressurization at constant injection rate (see ''Appendix'' A).
In this contribution, we assume that maximum over-pressure occurring at injection point x ¼ 0 remains always below the ambient effective normal stress r 0 o ¼ r n À p o applied on the fault plane, i.e. pð0; tÞ implying that the minimum principal effective stresses r 0 n ðx; tÞ remain always compressive (positive) throughout pressurization and hydraulic fracturing type of failure never occurs (which would require the full coupling between flow and elastic deformations). The along-fault pore-pressure diffusion changes the local effective normal stresses and hence drives the symmetric slip propagation when the Mohr-Coulomb criterion (3) is locally violated. Shear deformations, instead, do not affect pore-pressure evolution since shear-induced dilatancy or compaction is neglected and so fault permeability / porosity changes too. Although this assumption is debatable, since such inelastic deformations do affect slip stability on a planar fault with frictional weakening properties (Garagash and Rudnicki 2003;Zhang et al. 2005; Lecampion 2019), we want to minimise the complexities in the model and focus solely on the effect of injection rate ramp-up on the potential nucleation and arrest of dynamic slip. The hydromechanical model, therefore, is only one-way coupled and it is equivalent to the one proposed by Garagash and Germanovich (2012), with the difference that the point injection volumetric rate is not constant in time but changes according to (8).
We have introduced in the model an additional parameter t c , therefore we expect another dimensionless parameter governing the hydro-mechanical fault's response (on top of the ones introduced by Garagash and Germanovich (2012)). For sake of completeness we present and discuss in the next section all the dimensionless governing parameters resulting from scaling analysis (see ''Appendix'' B for more details), as well realistic values that are then used in the numerical simulations.

Dimensionless governing parameters
Upon normalization of all the governing equations (1-8) following Uenishi and Rice (2003); Garagash and Germanovich (2012) (see ''Appendix'' B), the hydromechanical fault response depends only on four dimensionless parameters: • Stress criticality s o s p , which represents the closeness of the ambient fault stress state to failure (and thus to the peak shear strength s p ¼ f p r 0 o ). Levandowski et al. (2018) claims that stress criticality is the most important factor for induced earthquakes hazard. Favourably oriented frictional weakening faults with respect to the in-situ stress field, typically characterised by a large stress criticality ( are very susceptible to host run-away seismic ruptures. Indeed, a little stress perturbation is sufficient to re-activate slip and its velocity propagation tends to diverge rapidly due to friction weakening and possibly other weakening mechanisms, such us flash-heating and thermal pressurization (Viesca and Garagash 2015). Garagash and Germanovich (2012) have shown via a stability analysis that critically stressed pressurized faults, for which the relation s o [ s r ¼ f r r 0 o is strictly satisfied, host always the nucleation of an unabated dynamic event. Such a run-away rupture, however, can be suppressed when shear-induced dilatancy kicks-off and dilatant hardening stabilises slip propagation (Lockner and Byerlee 1994;Segall et al. 2008;Ciardo and Lecampion 2019). Critically stressed faults in seismogenic zones have been observed in Oklahoma and Southern Kansas (Qin et al. 2019), in the German continental deep drillhole (KTB) (Ito and Zoback 2000) and in central California along the San Andreas fault system (Zoback et al. 1987;Rice 1992) to cite a few examples. On the other hand, fault zones not favourably oriented to the local in-situ stress field are characterised by low stress criticality ( s o s p J0), which implies that larger over-pressures are required to activate slip. By making the analogy with critically stressed pressurized faults exhibiting linear slipweakening behaviour, low stress criticality and s o \s r implies a quasi-static, stable slip propagation with eventually a nucleation and arrest of a dynamic event at large pressurization time (Garagash and Germanovich 2012). The role of the initial effective stress state on the nature of seismicity has been also investigated in fault lab experiments by Passelègue et al. (2020), where they show that faults with similar frictional properties can rupture at both slow and fast rupture velocity depending on their initial stress criticality, in agreement with theoretical predictions based on linear elastic fracture mechanics (LEFM).
• Friction weakening ratio the shear stress drop within the crack tips (with d r being the amount of slip at which the friction coefficient goes to zero if an unlimited linear slip-weakening friction law is considered). According to slip laboratory experiments of granite intact specimens, the slip weakening distance d r is approximately half a millimetre (Rice 1980), but it can drop even below 0.1 mm (Wong 1986).
Similarly, d w is on the order of fault's asperities and thus ranges between 0.1 and 10mm. The friction weakening ratio, therefore, can assume any values within the interval 0; 1 ð Þ.
Its value is typically in the order of a meter, but roughly one order of magnitude of variation is plausible due to the variation of d w with fault's roughness and r 0 o with hydrogeological conditions. In normally pressurized formations, typically located in the Earth's upper crust and characterised by lithostatic gradient and hydro-static pore-pressure conditions (Brace and Kohlstedt 1980; Grawinkel and Stockhert 1997), r 0 o may be a fraction of megapascals, while a boost to hundred of megapascals may be obtained in over-pressurized formations located below the fluid retention depth (Suppe 2014). Assuming a fluid viscosity parameter l $ 10 À3 Pa Á s (water), a fault gouge permeability k f in the order of $ 10 À16 m 2 (Wibberley et al. 2008) and hydraulic width w h in the order of meters (representative of a relatively-high permeable thick unit typically present in a damage zone of a fault zone-see Sect. 4.2 for a more indepth discussion), the injection rate characteristic scale Q Ã ranges between $ 10 À6 m 2 =s and 10 À3 m 2 =s. If we assume that the out-of-plane fault length is $ 1 km, then the volumetric injection rate scale would range between few to hundreds litres per second. The maximum injection rate Q m is a design parameter that can vary considerably from tens or hundreds litres per second in hydraulic-fracturing operations (see (Holland 2013) for one example where [ 160 l/s were injected in south-central Oklahoma), to few or fractions of litres per second for hydro-shearing stimulations (see for instance hydro-shearing experiments in Grimsel Test Site, Switzerland, where the maximum injection rate never exceeded 0.5 l/s (Villiger et al. 2020)).
The dimensionless parameter Q m Q Ã , therefore, can assume relatively low or large values depending on the specific problem configuration. Since hydraulic fracturing type of failure is not considered in our model, all the numerical results that will be presented later are characterised by a relatively • t c t w ! ratio between the ramp-up time scale and . This ratio may vary considerably due to the variation of slipping patch length-scale a w previously discussed, the specific fault hydraulic properties considered as well as the choice of the parameter t c . Assuming a compressibility parameter c f $ 0:5 GPa À1 (Wibberley 2002) and the hydraulic parameters previously defined (albeit the fault longitudinal permeability may range within the interval 10 À19 À 10 À15 m 2 ), we estimate an hydraulic diffusivity of a $ 2 Á 10 À4 m 2 =s, which leads to a diffusion time scale of t w $ 21 min for a w ¼ 1 m. The parameter t c may be in the order of several minutes for relatively fast fault pressurizations (see for instance one cycle of the injection protocol used for one hydro-shearing stimulation in Grimsel, Switzerland (Amann et al. 2018)), but it may get up to several hours or few days for relatively slow pressurizations (see for instance the injection protocol used for hydraulic stimulation of Pohang Enhanced Geothermal System, 2017 (Yeo et al. 2020;Hofmann et al. 2018)). The ratio t c =t w , therefore, may be very large, in the order 1 or possibly very small.
In the following we explore the parameter space identified by these dimensionless ratios via numerical simulations. We vary them systematically in order to investigate their impact on the nucleation and arrest (if occur) of dynamic fault slip.

Numerical results
We use a fast boundary element based solver in order to solve numerically the one-way coupled hydromechanical problem (1-11) (see details in (Ciardo et al. 2020)). It is suited for 2D non-linear geomechanical problems involving localized inelastic deformations along pre-existing structural discontinuities, such as faults or fractures. The elasto-static balance of momentum (1) is discretized using displacement discontinuity method that, together with the discretized form of shear weakening Mohr-Coulomb criterion (3), lead to a non-linear system of equations for the unknowns inelastic deformations (or displacement discontinuities) and effective tractions. For a given pore-pressure history calculated using (10) and (11) (with the latter evaluated numerically using a trapezoidal quadrature rule), such a resulting system is solved iteratively using fixed point iterations combined with under-relaxation (Quarteroni et al. 2000), and adopting a fully implicit integration scheme in time. A speed up of computations is obtained via hierarchical approximation of the fully populated elasticity matrix (on top of memory reduction), together with the use of an adequate block preconditioner that improves spectral properties of the resulting matrix of coefficients (see full details in (Ciardo et al. 2020)). All the numerical results that will be presented in the following are obtained by considering a sufficiently long planar fault such that L o =a w ¼ 40, with L o being its half-length. The planar fault is then uniformly discretized with 2 Á 10 3 equal-sized straight elements, resulting in having 50 elements inside the slip weakening region near crack tips (sufficient to accurately capture the model non-linearity).

Ultimately stable fault (s o \s r )
Firstly, we present the case of fluid injection into an ultimately stable fault characterised by s o =s r \1, where s r ¼ f r r 0 o is the residual frictional strength prior pressurization. This condition automatically implies that the fault is far from being critically stressed, and thus the ratio s o =s p is kept relatively low. By fixing the frictional weakening ratio to f r =f p ¼ 0:6, we investigate the fault's hydro-mechanical response for different values of normalized maximum injection rate Q m =Q Ã , in the two plausible limiting scenarios of t c =t w ( 1 and t c =t w ) 1.
3.1.1 Fast ramp-up: the limit when t c =t w ( 1 If the diffusion time-scale t w is much larger than the ramp-up time t c of injection rate, slip activation and propagation occurs when the injection rate has already reached its maximum value Q m =Q Ã . The quick rampup of injection rate has a negligible effect on porepressure evolution and the hydro-mechanical fault response is essentially driven by injection at constant volumetric rate. Garagash and Germanovich (2012) studied extensively the effect of such a type of fault pressurization on the nucleation and potential arrest of a dynamic event. They came up with a slip regimes  Garagash and Germanovich (2012)), as function of stress criticality s o =s p , friction weakening ratio f r =f p and normalized (constant) injection rate Q m =Q Ã . Region (1) corresponds to the scenario of an ultimately stable fault in which pressurization leads always to quasi-static (stable) slip propagation. Region (2) corresponds to the case of quasi-static slip propagation, followed by a nucleation and arrest of dynamic slip. Regions (4b) and (4a), instead, represent the scenarios of an unstable fault exhibiting quick quasi-static crack propagation followed by an unabated dynamic rupture, whose nucleation is affected or not by residual friction coefficient f r , respectively. Finally, in region (3) the nucleation of a run-away rupture on an unstable fault is preceded by a finite-sized seismic event diagram reported in Fig. 2, in which a frictional weakening planar fault may experience a finite-sized dynamic event (region (2)), a run-away rupture (affected or not by the residual friction coefficient f r , region (4b) and (4a) respectively), a finite-sized dynamic event followed by an unabated seismic rupture (region (3)), or only aseismic slip (region (1)). This depends on the particular set of dimensionless parameters that identifies a specific initial fault loading condition, pressurization and frictional property. Their theoretical and semi-analytical results, therefore, provide benchmark solutions for our numerical results.
Under ultimately stable conditions (i.e. s o \s r , left side of Fig. 2), Garagash and Germanovich (2012) proved via a stability analysis that a fault experiences a stable, quasi-static growth of slipping patch throughout sustained pressurization (i.e. aseismic slip). However, for sufficiently low injection rates Q m =Q Ã , the fault exhibits a finite-sized dynamic event, whose nucleation time t n is considerably larger than fluid diffusion time-scale t w (and thus t c ( t w ( t n ).
These considerations are also confirmed with our numerical results reported in Fig. 3, where the normalized time evolution of half crack length a=a w is displayed for two different values of Q m =Q Ã ¼ 0:05; 0:5 and different values of low/moderate stress criticality s o =s p (with the ratio t c =t w equal to 0.01 and s o always below or equal s r ).
For the case of low injection rate Q m =Q Ã ¼ 0:05, a nucleation of dynamic event followed by an arrest always occurs for each value of stress criticality considered, with larger dynamic run-out distances for increasing values of s o =s p (see Fig. 3-left). Since the ambient effective stress states applied on the fault plane are far from failure, resulting in a stable quasistatic slipping patch propagation before reaching the nucleation time t n (or in normalized form ffiffiffiffiffiffiffiffi ffi 4at n p =a w ), the shear crack tips are located well within the pressurized region. At nucleation time, the rupture is not affected by the residual friction coefficient f r (due to the low slip accumulation during quasi-static slip propagation phase) and its half-length assumes the asymptotic expression (Garagash and Germanovich 2012) A careful investigation of Fig. 4 that displays the corresponding slip, shear stress and friction coefficient profiles at different time snapshots (with dashed lines corresponding to the ones at nucleation time) and for two stress criticality values s o =s p ¼ 0:4; 0:6 confirms these theoretical predictions. It is worth noting at the Appendix'' C) associated with constant injection rate type of pressurization, which represents a good approximation of the hydro-mechanical fault response when t c =t w ( 1. In this particular example, indeed, the normalized ramp-up time is ffiffiffiffiffiffiffiffi 4at c p =a w ¼ 0:1 and thus shear crack activation and propagation is essentially driven at constant injection rate. Red dots, instead, denote the nucleation and arrest of dynamic slip arrest of the dynamic slip propagation, the friction coefficient has reached its residual value over almost the entire crack. The subsequent shear crack propagation, therefore, remains always stable in time as depicted in Fig. 3-left. In addition to this, Garagash and Germanovich (2012) proposed an asymptotic solution for the nucleation time t n that is valid for an ultimately stable fault that exhibits a dynamic event right after slip activation (similarly to the scenarios reported in Fig. 3-left). Under this condition, indeed, the shear crack at nucleation time is confined near injection point and is subjected to a uniform overpressure equal to the minimum value required to activate slip, i.e.
pðx ' 0; t n Þ ¼ ðs p À s o Þ=f p . This asymptotic solution reads (Garagash and Germanovich 2012) ffiffiffiffiffiffiffiffi ffi 4at n p where DP w ¼ Q m a w l 2k f ffiffiffi p p w h is the characteristic porepressure drop over the distance a w . In order to guarantee that crack instability follows shortly after slip activation, DP w must be very small compared to pðx ' 0; t n Þ, resulting in the following necessary condition The comparison between the nucleation times associated with numerical simulations reported in Fig. 3-left (for which Q m =Q Ã ¼ 0:05 strictly satisfies the condition (15) for each value of stress criticality considered) and the asymptotic solution (14) is displayed in Fig. 5. A good agreement is obtained, specially for larger values of s o =s p for which the assumption of nucleation after activation is truly valid (see Fig. 3-left). Obviously, lower values of Q m =Q Ã would lead to an earlier nucleation and thus a closer match between numerical and asymptotic solution.
For a larger value of injection rate Q m =Q Ã ¼ 0:5 (such to fall into region (1) of Fig. 2 for s o =s p \0:6), instead, the slipping patch propagates always quasistatically for all values of s o =s p , before reaching the maximum pressurization time ffiffiffiffiffiffiffiffiffiffiffiffi 4at max p =a w at which the condition (12) is violated (and at which all the while dashed vertical red lines denote the asymptotic solution (13) provided by Garagash and Germanovich (2012) simulations were stopped). One exception is a very small seismic event corresponding to the boundary case of s o =s p ¼ 0:6-see Fig. 3-right. For values of normalized time larger than ffiffiffiffiffiffiffiffiffiffiffiffi 4at max p =a w , a mixed mode fracture would propagate symmetrically from injection point, with a slipping patch located ahead of the tips of an opening (frictionless) patch. This problem has been studied by Azad et al. (2017) whose findings reveal that: i) the critical shear crack length a n =a w and slip at the hydraulic fracture tips increase rapidly with decreasing value of stress criticality s o =s p (with s o [ s r ), ii) the slip instability weakly depends on the injection rate applied and iii) fault slip development remains stable for s o s r . From these considerations we can infer that after ffiffiffiffiffiffiffiffiffiffiffiffi 4at max p =a w the propagation of the slipping patch ahead of the hydraulic fracture tips will always be stable, with an increasing yet quasi-static slip rate for values of stress criticality approaching the boundary value s o =s p ¼ 0:6.
These numerical results confirm that in case the maximum flow rate Q m is reached before pore fluid can actually diffuse into the fault, the hydro-mechanical fault response is essentially driven by pressurization at constant volumetric rate and a good match with theoretical results of Garagash and Germanovich (2012) is obtained. This is further strengthened by looking at the comparison between the numerical results reported in Fig. 3 (black lines) and the smallscale yielding (s.s.y.) asymptotic solutions (dashed grey lines) associated with pressurization at constant injection rate only (see ''Appendix'' C). For a=a w J2, indeed, a perfect match is obtained, suggesting that the effect of the quick ramp-up phase on slip propagation is negligible.
3.1.2 Slow ramp-up: the limit when t c =t w ) 1 We present here the opposite scenario in which injection rate increases linearly in time but the maximum value is reached after the pore fluid substantially diffuses along the fault plane (i.e. t c ) t w ). Under this condition, pore-pressure perturbation activates a shear crack (slip) during the ramp-up of injection rate, and the potential nucleation of a dynamic event would occur when pressurization is approaching the changing time t c . In other words, t c ) t w and t n .t c .
Scaling analysis reported in ''Appendix'' B suggests that the pressurization rate, i.e. how quick is the ramp-up of injection rate before time t c , does play a role in the hydro-mechanical fault response. In order to investigate its effect on the potential nucleation and arrest of dynamic slip on ultimately stable faults, we run several simulations keeping the ratio t c =t w constant (and much larger than 1) and varying the maximum injection rate Q m =Q Ã . This is equivalent to keeping the normalized injection rate constant and varying the (large) ratio t c =t w , resulting in relatively low or large pressurization rates. Fig. 6 displays the time evolution of half-crack length for different values of stress criticality, a ratio t c =t w ¼ 2 Á 10 3 and two maximum injection rates Q m =Q Ã ¼ 0:05; 1; representative of low and moderate fault pressurization rates, respectively. Notice that the maximum injection rates considered are similar to the ones used in the previous case of t c =t w ( 1 (see Fig. 3), although the pressurization rates are considerably different. From a comparison of Figs. 6-right and 3-right, we clearly observe that, in the case of t c =t w ) 1, a moderate maximum injection rate is not sufficient to quench the finite-sized dynamic slip event for each value of stress criticality (unlike the case previously discussed). This is certainly due to the different type of fault pressurization that drives the slipping patch expansion. In this case, indeed, slip is driven by linear increase branch of injection rate and pore-pressure at injection point evolves proportional to $ t 3=2 (see Eq. (10)). A further comparison of Figs. 3-left and 6-left for the exact same value of Q m =Q Ã ¼ 0:05 reveals that the dynamic run-out distances are rather similar for each corresponding value of stress criticality, while the nucleation times differ considerably. With the similar approach of Garagash and Germanovich (2012), we derive an asymptotic expression for the nucleation time that is valid when nucleation of dynamic slip follows shortly shear crack activation. By setting pðn ' 0; t n Þ from equation (10) equal to ðs p À s o Þ=f p , i.e. equal to the minimum over-pressure required to activate slip, we get the following asymptotic expression for the normalized nucleation time ffiffiffiffiffiffiffiffi ffi 4at n p where DP w is the same characteristic pore-pressure drop reported in Eq. (14). Unlike the asymptotic solution (14) valid for t c =t w ( 1, in this case the normalized nucleation time varies non-linearly with stress criticality s o =s p , with a power law exponent equal to 1/3. Furthermore, it appears the dependency on the dimensionless ratio t c =t w with the same power law exponent (as expected from scaling analysis), revealing that the nucleation time of dynamic slip does depend on how quick the injection rate ramp-up occurs in time, and thus on how large is the ratio  Fig. 7. A good agreement is obtained, with a closer match for larger values of s o =s p due to the earlier nucleations after slip activations (see Fig. 6-left). The asymptotic solution (16) gives also an estimation of the normalized position of the fluid front at nucleation time. A careful inspection of Fig. 6-left reveals that, at the onset of dynamic slip, the slipping patch front lags well within the pressurized region, and their relative distance decreases dramatically after the arrest of dynamic event (with fluid front always located ahead the slip front). Furthermore, we can observe that our numerical results (black lines) match  perfectly with the small-scale yielding (s.s.y.) asymptotic solutions associated with linear increase of injection rate (dashed grey lines-see ''Appendix'' C for more details) for each value of stress criticality. Since the extent of the arrested dynamic crack is always much larger than the slipping patch length scale a w (thus s.s.y. condition is truly valid) and the corresponding nucleation time is known analytically from (16), we can calculate analytically the dynamic run-out distances a arr directly from the small-scale yielding solution resolved at nucleation time. By simply replacing the nucleation time t n obtained from Eq. (16) into the increment of stress intensity factor associated with ramp-up of injection rate (43), we can solve the resulting implicit equation obtained from propagation criterion (46) for the unknown arrested crack lengths. Obviously, this analytical equation is only valid in the case of early dynamic crack nucleation after activation, but it allows to investigate systematically the effect of the other dimensionless parameters on the arrest of dynamic rupture.
We examine the effect of pressurization rate by fixing the maximum injection rate to Q m =Q Ã ¼ 0:05 (such to satisfy the necessary condition in Eq. (16) for the whole range of stress criticality s o =s p ) and varying the large ratio t c =t w , obtaining thus a relatively quick or slow ramp-up of injection rate. In Fig. 8, we display the normalized dynamic run-out distances a arr = ffiffiffiffiffiffiffiffi ffi 4at n p in function of stress criticality and friction weakening ratio f r =f p . We can observe that, for a given ramp-up of injection rate, the nucleation of a dynamic event is expected on a wider range of stress criticality for intermediate values of friction weakening ratio (compared to the cases of f r =f p ! 0 or f r =f p ! 1). Moreover, for a given value of f r =f p , the range of stress criticality in which a dynamic event is expected increases for increasing values of t c =t w (i.e. for decreasing values of Q m t w t c Q Ã ), revealing that low pressurization rates promote the nucleation of a finitesized dynamic rupture on ultimately stable faults. It is also interesting to note that the arrest of the dynamic Fig. 8 Normalized dynamic run-out distances a arr = ffiffiffiffiffiffiffiffi ffi 4at n p in function of stress criticality s o =s p and friction weakening ratio f r =f p , for different ramp-up scenarios of injection rate, i.e. different values of Q m t w t c Q Ã . These latter are obtained by fixing the maximum injection rate at Q m =Q Ã ¼ 0:05 and varying the (large) ratio t c =t w (notably t c =t w ¼ 2 Á 10 3 ; 2 Á 10 4 ; 2 Á 10 5 ) 1 ). In all the cases nucleation occurs prior reaching time t c , even for low values of f r =f p rupture always occurs within the pressurized region, i.e. a arr ffiffiffiffiffiffiffiffi ffi 4at n p .1, regardless of i) how quick the injection rate increase in time, ii) friction weakening ratio considered and iii) stress criticality value. Furthermore, for a given value of f r =f p and s o =s p , the dynamic run-out distance increases for decreasing values of Q m t w t c Q Ã , i.e. for slower ramp-up of injection rate. This can be grasped clearly from Fig. 9 where the normalized ramp-up rate Q m t w t c Q Ã is plotted against the normalized dynamic run-out distance a arr =a w in a log-log plot. As one can see, the arrest of the dynamic slipping patch decreases non-linearly with increasing values of ramp-up rates, approximately proportional to an inverse power-law with exponent equal to -0.472 (see red line in Fig. 9).

Unstable fault (s o [ s r )
Finally, we present the case of injection into an unstable fault characterised by s o =s r [ 1 at ambient conditions. In this particular scenario the fault is critically stressed and thus prompt to fail. A small strength perturbation, due to for instance a small porepressure increment p, always leads to a quick quasistatic shear crack propagation followed by a run-away dynamic rupture (albeit in some circumstances a finite-sized dynamic event may precede the run-away rupture-see region (3) of Fig. 2 valid for the injection scenario at constant volumetric rate). Assuming a shear crack activation during the early ramp-up of injection rate (valid when t c is not much smaller than t w and thus pore fluid can diffuse within the fault plane), the slipping patch outpaces rapidly the pore fluid front and, at nucleation time t n , the following condition hold a n a w ) with t n much smaller than both fluid diffusion timescale t w and ramp-up time scale t c . The pressurized region at nucleation time, therefore, is always confined near injection point and the corresponding porepressure distribution can be approximated using an equivalent point force distribution pðx; tÞ ' DPðtÞd dirac ðxÞ, with DPðtÞ defined in (10).
Based on the previous work of Uenishi and Rice (2003), Garagash and Germanovich (2012) showed that there exist an asymptotic solution in terms of critical shear crack length that is universal (i.e. independent of the particular type of injection scenario) and it reads a n a w ' 0:579 In addition to this, they developed semi-analytically an outer and inner asymptotic solutions that are valid at different fault position with respect to fluid front location (see ''Appendix'' D). The outer solution is so called because is valid outside the pressurization region, i.e. for x j j ) ffiffiffiffiffiffi ffi 4at p . The inner asymptotic solution, instead, is valid for x j j. ffiffiffiffiffiffi ffi 4at p , i.e. near injection point due to the limited extension of pressurization region on unstable faults before instability.
Since the outer asymptotic solution is universal, i.e. it is valid for any type of peak pore-pressure distribution (see ''Appendix'' D), we can make use of such a solution to derive an asymptotic expression for the nucleation time t n in function of the ramp-up of injection rate (and compare it with the solution of Garagash and Germanovich (2012) valid for pressurization at constant flow rate that can actually be retrieved here in the limit of t c ! 0). Based on the outer solution, indeed, the integrated net over-pressure along the fault at crack instability is approximately given by (Garagash and Germanovich 2012) Fig. 9 Normalized dynamic run-out distances a arr =a w (black dots) in function of pressurization rate Q m t w t c Q Ã (with t c =t w ) 1), for an ultimately stable fault characterised by s o =s p ¼ 0:55 and friction weakening ratio of f r =f p ¼ 0:6. The red solid line represents a linear fit of the numerical results (black dots) in the log-log plot where P ' 0:8369 is the scaled magnitude of the point force, independent of the type of fault pressurization. By replacing Eq. (10) into Eq. (19), after some algebra, we obtain ffiffiffiffiffiffiffiffi ffi 4at n p where DP w is the same characteristic pore-pressure drop of Eq. (14). It is interesting to note that the normalized nucleation time in the case of an unstable fault subjected to a ramp-up of injection rate varies non-linearly with both stress criticality s o =s p and pressurization rate Q m Á t w Q Ã Á t c (with a power law exponent equal to 1/4). A comparison with the normalized nucleation time associated with constant injection rate type of pressurization (Garagash and Germanovich 2012) suggests that the ramp-up of injection rate prior the maximum plateau may promote an earlier or later nucleation of a run-away dynamic event, with respect to the injection scenario at constant volumetric rate. Indeed, by taking the ratio between equation (20) and equation (21), we obtain 1:20636 Á f p DP w s p À s o t c t w 1=4 , implying that an earlier nucleation is expected when for a given maximum value of injection rate Q m .
In order to verify all these theoretical predictions, we run several numerical simulations with different stress criticality values, such to obtain highly critically stressed fault conditions for which condition (17) is strictly satisfied. We fixed the maximum injection rate at Q m =Q Ã ¼ 100, the friction weakening ratio at f r =f p ¼ 0:6 and the normalized changing time at t c =t w ¼ 0:5. As we can observe from Fig. 10, the quick quasi-static shear crack propagation is always followed by a run-away dynamic rupture, whose nucleation time is in good agreement with theoretical prediction of Eq. (20) (with a better accuracy for larger values of s o =s p due to a more strict validity of condition (17)). At the instability, the normalized crack length a=a w is approximately $ 0:579 for each value of stress criticality considered (as expected) and the corresponding scaled slip distributions match very well the outer and inner asymptotic solutions at x j j ) ffiffiffiffiffiffi ffi 4at p and x j j. ffiffiffiffiffiffi ffi 4at p , respectively (see Figs. 12 and 13 in ''Appendix'' D).

Discussions
In this section, we first discuss the limitations of our modelling approach and their effects on nucleation and arrest of dynamic fault slip. We then talk about the fault zone permeability structure and its relation to the parametrization of the model. Finally, we present the implications of our results on injection-induced seismicity, a topic that has raised the attention of scientific community and public opinion since second half of 20 th century (Simpson 1986).

Model limitations
The most severe approximation in our model is the neglect of dilatancy during fluid-driven shear crack propagation. Strong dilatant behavoir associated with sliding over fault's asperities has been observed in laboratory experiments (Lockner and Byerlee 1994;Samuelson et al. 2009) as well as during field experiments in the context of geothermal energy exploitation (Batchelor 1985). Under the assumption of no fluid leak-off in the relative undamaged rock around the fault, shear-induced dilatancy does affect pore-pressure evolution, which in turn leads to a feedback on slip propagation (due to full coupling between flow and elastic deformations). Under undrained conditions, indeed, shear-induced dilatancy leads to a pore-pressure drop and thus to a local increase of effective normal stress (dilatant hardening) (Rudnicki 1979;Segall and Rice 1995). Ciardo and Lecampion (2019) have shown that such a dilatant hardening effect, quantified by a scaled undrained pore-pressure drop Dw h w h c f r 0 o with Dw h being the increment of fault opening, does impact the transition between aseismic and seismic slip on frictional weakening immature faults. They showed, in fact, that a large dilatant behaviour suppresses the nucleation of run-away seismic rupture on otherwise unstable faults, even for sustained increases of fault permeability with slip. On ultimately stable faults, instead, dilatant hardening effect delays the nucleation of finite-sized seismic event and increases its dynamic run-out distance. We can certainly say, however, that the findings obtained in this work are valid for a fault whose dilatant compliance is relatively small, i.e. for In this contribution we also neglect other weakening mechanisms that can kick in during the onset of dynamic crack propagation, such as thermal pressurization and flash heating (Rice 2006).
Thermal pressurization of pore-fluid by rapid shear heating of fault gouge unit has been showed to be a prominent process of fault weakening (Viesca and Garagash 2015). This mechanism depends on fluid thermodynamics properties and drives the shear strength loss at low fluid pressure conditions (Acosta et al. 2018). Garagash and Germanovich (2012) showed that dynamic weakening due to thermal pressurization increases the dynamic run-out slip distances on ultimately stable faults. However, they also stated that such an effect is relatively small, due to the fact the most of dynamic weakening is expected to occur at slip scale d w;dyna that is likely larger than d w . Without loss of generality, we can claim that the results presented in this contribution in terms of dynamic run-out distances on ultimately stable faults are valid for sufficiently low values d w compared to the characteristic dynamic slip-weakening distance d w;dyna (which depends on heat properties of both fault gouge unit and injected fluid).
Flash heating on fault's asperities, instead, kicks in when the slip velocity exceeds $ 0:1m/s, a scenario that is very plausible during dynamic crack propagation, which characteristic slip rate at the rupture tip typically exceeds $ 10m/s (Garagash 2011). Laboratory experiments on various types of rocks have shown a drop of friction coefficient of up to one order of magnitude due to flash heating on asperities contacts (Di Toro et al. 2011). These considerations suggest that this weakening mechanism would impact the dynamic run-out distances presented in this contribution. However, this is outside the scope of this work and it is left for future investigations.
Finally, we adopted a simple linear slip-weakening friction law compared to a more elaborate Rate-and State-friction model (Dieterich 1979). Although a number of studies have demonstrated that the linear slip-weakening friction model is a good approximation of the velocity weakening R&S friction law (Uenishi and Rice 2003;Rubin and Ampuero 2005;Garagash 2021) even at slip instability (Viesca 2016b, a), it does not allow to investigate scenarios in which fault frictional properties evolve during pressurization, for instance from velocity weakening to velocity strengthening depending on current slip velocity (as observed in laboratory experiments by Cappa et al. (2019)). Indeed, this can only be captured by using R&S friction law with heterogeneous distribution of a and b parameters. Future works with incorporation of R&S friction model will follow.

Fault zone permeability structure
Fault zones are complex lithological zones that typically contain a fault core with possibly branching subsidiary faults and a damage zone (Faulkner et al. 2010). Their geological structures are highly heterogeneous, anisotropic and discontinuous, and thus characterised by highly variable permeability structures (Lockner et al. 2009). A common description of a fault zone is a central core surrounded by a damage zone. Exhumed fault zones typically reveal the presence of several relatively higher permeability units compared to the host rock, with values around k f $ 10 À16 À 10 À15 m 2 (Wibberley 2002)). In case of an injection in or near the damage zone, the pore pressure would distribute in a relatively large region, as the damage zone itself would plausibly serve as main fluid conduit. Generally, the permeable units in the damage zone (fractures/faults) are relatively thick, in the order of meters, and cover a wide range of length scales. On the other hand, the fault core is typically thinner (centimetres/millimetres) and composed of gouge, cataclasite and/or ultracataclasite. Its permeability is relatively low, in the range of k f $ 10 À19 À 10 À18 m 2 (Wibberley 2002), and therefore it acts as an impermeable barrier for pore-fluid diffusion, at least at the time scale of injection. However, the core is the part of the fault zone that accommodates most of the slip (Faulkner et al. 2010).
The model presented in this contribution is a simplification of such a complex fault zone structure. We assume that fluid is injected into or near a highly permeable thick unit that accommodates both pore fluid diffusion and slip propagation. This assumption has an important consequence on the hydro-mechanical response since the permeability of the conductive unit k f and its thickness w h enter in the fluid diffusion time scale t w (only the former) as well as in the characteristic injection rate scale Q Ã . These two scales, indeed, are used to normalised the two operational parameters in the model, i.e. the ramp-up time scale t c and the maximum injection rate Q m , respectively.

Implications on injection-induced seismicity
The results above may indicate that injection-induced seismicity can be mitigated by controlling operational parameters, in particular the injection volumetric rate or, more specifically to the problem addressed in this contribution, the pressurization rate Q m t w t c Q Ã . Upon fluid injection into a specific fracture zone, possibly indicating of a well-developed (mature) fault zone, pore-pressure perturbation is likely to activate slip on favourably oriented (and critically stressed) conductive fractures within the damage zone. The initial stable slip propagation could be quickly followed by run-away seismic ruptures that in turn could trigger events on other fractures and propagate over the entire fault zone. Our results show that the onset of the contained micro-seismicity, however, may be delayed in time for slow increases of injection rate prior reaching a steady-state phase. If a specific stable principal fault plane is targeted for fluid injection, instead, a larger pore-fluid perturbation is required to activate slip. A ramp-up of the injection rate, in this case, strongly affects the slip propagation, in particular when the ramp-up time scale is much larger than the along-fault diffusion time scale. Although counter-intuitive, a fast ramp-up of injection rate on stable fault planes would reduce the possibility of triggering a larger finite-sized seismic event (with respect to the one that could be triggered if fault pressurization occurs at constant injection rate), which can potentially turn into an unabated rupture due to the activation of other dynamic weakening mechanisms (as discussed in Sect. 4.1).
The results are essentially in line with previous laboratory and numerical observations. Indeed, previous results highlighted how the initial state of effective stress (Gischig 2015;Passelègue et al. 2018) and essentially its relationship to the frictional behaviour (Larochelle et al. 2021) control whether a fault slip is confined to the pressurized region or runaways.
Similarly, the numerical analysis by Alghannam and Juanes (2020) reports how the pressurization rate may influence the seismic reactivation of a fault zone. Here we generalize both approaches by showing that the ruptures and its final behaviour (confined or runaway) are linked to both the initial state of stress as well as the initial ramp-up of the injection, which is then closely linked to the pressurization rate.
The effects of the controlling operational parameters on risk of induced-seismicity associated with the specific physic-based model presented in this contribution are summarized in Fig. 11. It displays the effects of the initial ramp-up of injection rate on the different slip regimes (the same of Fig. 2), as function of all the other governing dimensionless parameters. From the characteristics of fault zone permeability structures discussed in the previous section, we deduce that the characteristic injection rate scale Q Ã might be realistically in the order of tens or even hundreds of litres per seconds due to the presence of highly permeable thick units that serve as fluid conduits. This leads to values of Q m =Q Ã relatively low for hydroshearing type of stimulations, implying that the risk of triggering a seismic event is also plausible in the case of injection into ultimately stable fault zones. This scenario is even more plausible for very large values of normalized ramp-up time scale t c =t w , which mean lower pressurization rates Q m t w t c Q Ã and thus a larger region (2) in Fig. 11. For slightly unstable faults (s o Js r ), instead, a lower pressurization rate would lead to a larger slip accumulation during the quick quasi-static phase of crack propagation. The friction coefficient, therefore, would drop quicker to its residual value and the probability of triggering a runaway dynamic rupture preceded by a finite-sized event is higher, resulting in a larger region (3). This also implies that the nucleation time of such a run-away rupture increases for lower pressurization rates, as more and more part of the shear crack has reached the residual shear strength during its quasi-static stable propagation.

Conclusions
In this study, we have extended the model of Garagash and Germanovich (2012) to account for an initial ramp-up of injection rate in time before the maximum plateau and investigated its effect on nucleation and arrest of dynamic fault slip. Despite the simplicity of the homogeneous model (planar bidimensional fault, uniform stress conditions and rock properties, linear weakening of friction, no dilatancy), it allows to get insight into the mechanisms that govern the transition between aseismic and seismic slip on a pressurized planar fault.
We have approximated standard injection protocols that consist of an initial step-wise ramp-up of injection rate in time with a linear increasing function, followed by a maximum plateau after a given time t c . We have solved numerically the coupled hydro-mechanical problem and explored the different slip regimes identified via scaling analysis. Our results show that the initial ramp-up of injection rate affects the slip Fig. 11 Map of slip regimes for a frictional weakening pressurized fault in which injection rate increases linearly in time, up to reach a maximum plateau after time t c . This parametric diagram is function of stress criticality s o =s p , friction weakening ratio f r =f p , normalized maximum injection rate Q m =Q Ã and normalized ramp-up time scale t c =t w . For a given value of Q m =Q Ã , the larger is the ratio t c =t w , i.e. the lower is the pressurization rate Q m t w t c Q Ã , the more wide are region (2) and (3). All the slip regimes of this Figure are the same of those reported and described in Fig. 2 propagation on ultimately stable faults (s o \s r ) if and only if the ramp-up time scale t c is much larger than fluid diffusion time scale t w . Notably, we have shown that slip stability is governed by the pressurization rate (quantified by the dimensionless parameter Q m t w Q Ã t c ) and thus by how quickly the injection rate increases before reaching time t c . From our results we can conclude that low pressurization rates applied on ultimately stable faults with frictional weakening properties and t c ) t w promote the nucleation of a finite-sized dynamic event. Moreover, the lower is the pressurization rate, the larger is the dynamic runt-out slip distance and hence the larger is the magnitude of the induced seismic event. We have also developed an asymptotic solution in terms of nucleation time that is valid when slip instability follows shortly crack activation and verified it with numerical simulations. When t c ( t w , instead, the initial ramp-up of injection rate can be neglected on ultimately stable faults and slip are driven by pressurization at constant volumetric rate. A good agreement with theoretical predictions of Garagash and Germanovich (2012) valid for that particular injection condition has been obtained.
Finally, we have demonstrated that the ramp-up of injection rate does affect the quick slip propagation on critically stressed faults prior the nucleation of a runaway dynamic rupture. By using the universal outer asymptotic solution of Garagash and Germanovich (2012), we developed a new analytical solution for the nucleation time and verified it with numerical simulations. This solution reveals that the initial ramp-up of injection rate may lead to an earlier or later nucleation of unabated dynamic event, compared to the case of fault pressurization at constant volumetric rate.
Acknowledgements This work has been mainly subsidised through the ERANET Cofund GEOTHERMICA (Project No. 200320-4001) by the Swiss Federal Office of Energy (SFOE), which is supported by the European Union's HORIZON 2020 programme for research, technological development and demonstration. The work was also supported by the ERC Synergy grant FEAR (Grant agreement No. 856559). Technical review comments by Elias Heimisson at SED are greatly appreciated. The authors acknowledge review comments from an anonymous reviewer. The analytical formulae and numerical methods described in the main text and appendices are sufficient to reproduce all the results presented in the paper.

Author
Contributions F.C. contributed to the Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Visualization, Writing original draft. A.P.R. contributed to Editing, Supervision.
Funding Open access funding provided by Swiss Federal Institute of Technology Zurich.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Appendix A: Along-fault pore-pressure diffusionanalytical solution Equation (5) or equivalently equation (9) is a parabolic, linear second-order partial differential equation with constant coefficients. It governs the one-dimensional over-pressure diffusion pðx; tÞ ¼ pðx; tÞ À p o inside the fault gouge unit with constant hydraulic diffusivity a L 2 =T ½ . Together with the boundary conditions pðAE1;tÞ ¼ 0; and the specific initial conditions, it represents a wellposed diffusion problem that can be solved analytically (Carslaw and Jaeger 1959).
Firstly, we solve the diffusion problem for the case of linear ramp-up of injection rate valid for t t c , for which the initial condition reads pðx; 0Þ ¼ 0 We use Laplace transform in time (with Laplace parameter s) in order to turn equation (9) into an ordinary differential equation (with transformed variables denoted with an hat^). With the initial condition (24), such a subsidiary equation reads o 2p ðx; sÞ ox 2 À k 2p ðx; sÞ ¼ 0 with whose analytical solution iŝ pðx; sÞ ¼ c 1 ðsÞ Á e kx þ c 2 ðsÞ Á e Àkx ð26Þ c 1 ðsÞ and c 2 ðsÞ are two constants that can be obtained in closed form using the boundary conditions (23) associated only with linear ramp-up of injection rate: characteristic scales in order to normalize elasticity equation (1) and shear stress evolution within the crack tips (3): sðx; tÞ ¼ s p Á T x; t ð Þ; dðx; tÞ ¼ d w Á D x; t ð Þ; pðx; tÞ ¼ r where d w is the slip weakening length-scale (see Fig. 1), s p ¼ f p r 0 o is the peak shear strength at ambient conditions and a is half-length of the shear crack. The corresponding governing equations upon introduction of (33) read: a w a slipping patch is driven essentially by pressurization at constant injection rate (for instance in the case of an ultimately stable fault in the limit when t c =t w ( 1), equation (46) is solved using the stress-intensity factor obtained in (44). Instead, when the shear crack is driven up to the nucleation of a dynamic event by the ramp-up of injection rate (for instance in the case of an ultimately stable fault with t c =t w ) 1), then equation (46) is solved using (43) (obviously up to the maximum value t c ). The comparison between the numerical results (black lines) and small-scale asymptotic solutions (dashed grey lines) in Figs. 3 and 6 shows a good match for slipping patches a larger than $ 2a w .
In order to check the accuracy of these theoretical predictions, we display in Fig. 12 the corresponding slip profiles at nucleation time t n of the numerical simulations reported in Fig. 10. As we can observe, the scaled slip distributions d n =ðd w Þ, which collapse into nearly one black line due to the scaling adopted, match the outer asymptotic solution (50) (denoted by a red solid line) for x=a n j j) ffiffiffiffiffiffiffiffi ffi 4at n p =a n , i.e. outside the pressurization regions whose extents are identified by blue dashed lines (for a given value of stress criticality s o =s p ). For x=a n ( 1, instead, the scaled slip distributions tend to converge to the near-field asymptotic expansion of the outer solution (52) (as expected).

Inner solution
Within the small pressurization region, i.e. for x.
ffiffiffiffiffiffi ffi 4at p , the normalized pore-pressure distribution is approximately pðxÞ=r 0 o $ 1, which means that the stress perturbation can be written as (Garagash and Germanovich 2012) for ( 1. Using this condition, Garagash and Germanovich (2012) showed that the normalized slip distribution at instability d n =d w for a particular porepressure profile WðsÞ is given as where d n ð0Þ is the critical slip at n ¼ x 4at ¼ 0, which can be obtained by matching the outer (x ) ffiffiffiffiffiffi ffi 4at p ) and inner (x. ffiffiffiffiffiffi ffi 4at p ) asymptotic solutions at intermediate distances (Garagash and Germanovich 2012) where ¼ ðs p À s o Þ=f p D pðt n Þ and the constant C is defined as Unlike the outer asymptotic solution, the inner solution does depend on the particular pore-pressure profile WðsÞ. By replacing the instantaneous spatial distribution (10) associated with ramp-up of injection rate, we obtain C ¼ 3:75576, allowing us to calculate the normalized slip distribution using (55) and (54) (with numerical evaluation of the integral R 1 À1 WðsÞln 1 À n=s j j ds). For the same injection condition, friction weakening ratio and initial loading conditions of the simulations reported in Fig. 10, we show in Fig. 13 the corresponding normalized slip distributions d n =d w at nucleation time t n on a linearlog plot. We can observe that, for each value of (large) stress criticality s o =s p , the corresponding inner Fig. 12 Corresponding normalized slip distributions d n =ðd w Þ (black solid lines) at nucleation time t n of the numerical simulations reported in Fig. 10. The red solid line denotes the outer asymptotic solution valid outside the pressurization regions (whose extents at instability are identified by blue dashed lines), i.e. for x a n ) ffiffiffiffiffiffiffiffi ffi 4at n p a n , while the grey solid line represents the near-field expansion of the outer asymptotic solution that is valid for x=a n ( 1 Fig. 13 Linear-log plot showing the corresponding normalized slip distributions d n =d w (black solid lines) at nucleation time t n of the numerical simulations reported in Fig. 10. The green solid lines denote the inner asymptotic solutions valid within the pressurization regions, i.e. for x ffiffiffiffiffiffiffiffi ffi 4at n p .1 asymptotic solution matches well the numerical results (back solid lines) for x j j= ffiffiffiffiffiffiffiffi ffi 4at n p .1 (as expected).