Three-dimensional fluid-driven stable frictional ruptures

We investigate the quasi-static growth of a fluid-driven frictional shear crack that propagates in mixed mode (II+III) on a planar fault interface that separates two identical half-spaces of a threedimensional solid. The fault interface is characterized by a shear strength equal to the product of a constant friction coefficient and the local effective normal stress. Fluid is injected into the fault interface and two different injection scenarios are considered: injection at constant volume rate and injection at constant pressure. We derive analytical solutions for circular ruptures which occur in the limit of a Poisson’s ratio ν = 0 and solve numerically for the more general case in which the rupture shape is unknown (ν 6= 0). For an injection at constant volume rate, the fault slip growth is self-similar. The rupture radius (ν = 0) expands as R(t) = λL(t), where L(t) is the nominal position of the fluid pressure front and λ is an amplification factor that is a known function of a unique dimensionless parameter T . The latter is defined as the ratio between the distance to failure under ambient conditions and the strength of the injection. Whenever λ > 1, the rupture front outpaces the fluid pressure front. For ν 6= 0, the rupture shape is quasi-elliptical. The aspect ratio is upper and lower bounded by 1/(1 − ν) and (3 − ν)/(3 − 2ν), for the limiting cases of critically stressed faults (λ 1, T 1) and marginally pressurized faults (λ 1, T 1), respectively. Moreover, the evolution of the rupture area is independent of the Poisson’s ratio and grows simply as Ar(t) = 4παλt, where α is the fault hydraulic diffusivity. For injection at constant pressure, the fault slip growth is not self-similar: the rupture front evolves at large times as ∝ (αt)(1−T )/2 with T between 0 and 1. The frictional rupture moves at most diffusively (∝ √ αt) when the fault is critically stressed, but in general propagates slower than the fluid pressure front. Yet in some conditions, the rupture front outpaces the fluid pressure front. The latter will eventually catch the former if injection is sustained for a sufficient time. Our findings provide a basic understanding on how stable (aseismic) ruptures propagate in response to fluid injection in 3-D. Notably, since aseismic ruptures driven by injection at constant rate expands proportionally to the squared root of time, seismicity clouds that are commonly interpreted to be controlled by the direct effect of fluid pressure increase might be controlled by the stress transfer of a propagating aseismic rupture instead. We also demonstrate that the aseismic moment M0 scales to the injected fluid volume V as M0 ∝ V .


Introduction
Fluid-driven frictional ruptures play an important role in earthquake and fault mechanics and can occur either as a natural process or be induced by human activities. Some examples of the natural source are related to fault valving behavior (Sibson 1992, Zhu et al. 2020) and metamorphic dehydration reactions (Wong et al. 1997, Hacker et al. 2003, Kato et al. 2010 in fault systems, whereas seismic swarms (Parotidis et al. 2005, Chen et al. 2012, Ross et al. 2020, Hatch et al. 2020) and aftershock sequences (Bosl & Nur 2002, Miller et al. 2004, Hainzl et al. 2016, Ross et al. 2017, Miller 2020) are other natural phenomena commonly attributed to the migration of fluids in fault zones. On the other hand, anthropogenic fluid injections associated with hydrocarbon and geothermal operations routinely produce microseismicity and have been extensively linked to the reactivation of faults (Healy et al. 1968, Deichmann & Giardini 2009, Keranen et al. 2013, Eyre et al. 2019, Ellsworth et al. 2019.
Evidence for fluid-driven frictional ruptures is generally inferred from the observation of seismicity spreading away from natural or human-related fluid injections in the Earth's crust (Shapiro et al. 1997, Parotidis et al. 2005, Hainzl et al. 2016, Ross et al. 2017, Goebel & Brodsky 2018, Ross et al. 2020. Observed seismicity is the result of unstable (seismic) frictional sliding that radiates detectable seismic waves. Nevertheless, seismic slip is not the only possible result of fluid injection. In fact, stable (aseismic) slip, which is more difficult to detect due to the virtual absence of elastodynamic waves, is a likely frequent result of the injection of fluids as demostrated by past large-scale fluid injections in the field (Hamilton & Meehan 1971, Scotti & Cornet 1994, Bourouis & Bernard 2007, recent laboratory and in-situ experiments (Guglielmi et al. 2015, Scuderi & Collettini 2016, Cappa et al. 2019, Passelègue et al. 2020, and recent cases of injection-induced seismicity (Wei et al. 2015, Chen et al. 2017, Eyre et al. 2019. As suggested by a number of recent experimental and observational studies (Wei et al. 2015, Guglielmi et al. 2015, Duboeuf et al. 2017, Eyre et al. 2019, Cappa et al. 2019, Bhattacharya & Viesca 2019, injection-induced aseismic slip may trigger seismicity by the transfer of solid stresses to unstable patches in pre-existing structural discontinuities, such as fractures and faults. Such spatio-temporal perturbation of the stress field is due to a quasi-statically expanding fluid-driven slipping patch that propagates along an initially locked and predominantly frictionally-stable pre-existing discontinuity. In this view and in the framework of this model, seismicity can be conceptually understood as the result of instabilities triggered by perturbating the pre-injection stress state of frictionally-unstable patches present either in the same pre-existing discontinuity (due to heterogeneities in rock frictional properties) or in others nearby the propagating rupture. The potential prominence of this triggering mechanism has increased in recent times since new investigations have suggested that fluid-induced aseismic slip can outpace the diffusion of fluid pressure (Eyre et al. 2019, Bhattacharya & Viesca 2019 and may be in fact the primary cause of observed seismicity during in-situ experiments (Guglielmi et al. 2015, Duboeuf et al. 2017) and responsible for the triggering of hydraulic fracturing-induced earthquakes (Eyre et al. 2019).
Recent efforts for understanding injection-induced aseismic slip have been motivated mostly by the sudden increase of seismicity due to anthropogenic fluid injection (Keranen et al. 2014, Bao & Eaton 2016, Goebel & Brodsky 2018. Nevertheless, understanding the mechanics of fluid-driven aseismic slip is indeed relevant to any phenomenon that is predominantly characterized by stable frictional slidding and the pressurization of interfacial fluids. This might be the case of, for instance, some seismic swarms (Chen et al. 2012, Hatch et al. 2020, aftershock sequences (Ross et al. 2017), and slow slip transients near the base of the seismogenic zone due to fault valving (Zhu et al. 2020) or metamorphic dehydration reactions (Kato et al. 2010).
Despite the apparent relevance of fluid-driven aseismic ruptures in a wide variety of natural and anthropogenic phenomena, the spatio-temporal evolution of aseismic slip in response to fluid injection remains poorly constrained in 3-D. This is, in part, due to the challenge of solving such a moving boundary value problem in which both fault slip and rupture shape are unknown. In this article, we investigate the mechanics of injection-induced fault slip by solving the problem of a fluid-driven frictional shear crack that propagates in mixed mode (II+III) on a planar fault that separates two identical half-spaces of a three-dimensional, linear elastic, and impermeable solid. The fault interface is saturated by pressurized fluid and it is characterized by a constant hydraulic transmissivity and a shear strength that is determined by the product of a constant friction coefficient and the local effective normal stress. We consider that fluid is injected into the fault interface under two injection scenarios: injection at constant volume rate and injection at constant pressure. The model is an extension to 3-D of a previous 2-D model presented by Bhattacharya & Viesca (2019). In the process, we also investigate the fundamental problem of crack-shape selection of a frictional shear crack under localized (point-force-like) and distributed effective shear loadings, including its dependence on the Poisson's ratio of the bulk. This paper is organized as follows. In section 2, we present the mathematical formulation of the problem and the chosen numerical methods. In section 3, we solve the problem of a stable rupture driven by injection at constant volume rate, for which we first derive an exact analytical solution in the case of axisymmetric circular ruptures, and then solve numerically for the more general case in which the x z f a u l t p l a n e f a u l t p l a n e rupture front, nominal pore pressure front, L(t) Q(t) r θ y Figure 1: Model schematic. A planar fault separates two semi-infinite linear elastic, homogeneous, and isotropic solids. The fault is characterized by a constant friction coefficient and is embedded in uniform initial fluid pressure and stress fields. Fluid is injected right into the fault through a wellbore located along the z axis. Fluid flow is confined within the fault plane and is fault-parallel and axisymmetric with regard to the z-axis. A quasi-static rupture of front R(t) whose shape is unknown is driven by axisymmetric fluid pressure diffusion that is characterized by a nominal fluid pressure front L(t).
rupture shape is part of the solution. Following a similar approach, we solve in section 4 the problem of a stable rupture driven by injection at constant pressure. Finally, section 5 discusses our findings and its possible implications to a number of fluid-driven earthquake-related phenomena such as injection-induced seismicity, seismic swarms, aftershock sequences and slow slip events.

Problem formulation
We consider a fault plane Γ located along z = 0 that separates two semi-infinite, homogeneous, isotropic and linear elastic solids (see Fig. 1). The fault interface is governed by Coulomb's friction with a constant friction coefficient. The initial stress tensor is uniform and is characterized by a shear stress τ o resolved on the fault plane that acts along the x direction, and a total normal stress σ o to the fault plane (that acts along the z direction). We assume that fluid is injected into the fault plane through, for instance, a wellbore, that is located along the z axis. We also assume that the solid is impermeable and the fault interface has a uniform and constant hydraulic transmissivity; fluid flow thus occurs only within the fault plane and is fault-parallel and axisymmetric with regard to the z axis. Owing to the planarity of the fault and the uniform direction of the initial shear stress τ o , fluid flow induces fault slip δ and changes in the shear stress τ resolved on the fault plane that are both characterized by a uniform direction along the x axis. The magnitude of the fault slip δ is maximum at the origin (the injection point) and vanishes along the rupture front R(t) = {(x, y) : δ(x, y, t) = 0}. The rupture front R(t) is unknown a priori and is to be determined as part of the solution.
The quasi-static elastic equilibrium that relates the fault slip δ to the shear stress τ on the fault plane Γ can be written as the following boundary integral equation (Hills et al. 1996) τ (x, y, t) = τ o + Γ K(x − ξ, y − ζ; µ, ν)δ(ξ, ζ, t)dξdζ, where τ o is the initial shear stress, µ is the shear modulus of the solid, ν is the Poisson's ratio, and K is the hypersingular (of order 1/r 3 ) elastostatic traction kernel (Hills et al. 1996). We adopt the convention of slip positive in clockwise rotation, δ(x, y, t) = u x (x, y, z = 0 where u x is the displacement component in the x direction. We also adopt the convention of normal stress positive in compression.
The fault is assumed to obey a Mohr-Coulomb shear failure criterion without any cohesion: where f is the constant friction coefficient, σ o − p(r, t) denotes Terzaghi's effective normal stress to the fault plane with p(r, t) being the spatiotemporally evolving fluid pressure, which is axisymmetric about the point of injection {O; r, θ, z} (see Fig. 1). We further assume that fault slip occurs without any normal displacement discontinuity, neither dilatant nor contractant, and that fault slip does not impact fluid flow.
We make the assumption that the surrounding rock can be considered impermeable compared to the fault itself at the scale of the injection duration. Single phase porous media flow (Bear 2013) in the fault thus reduces (after width averaging across the fault thickness w h ) to the following two-dimensional pressure-diffusion equation on the fault plane where S is a storage coefficient combining the effect of both fluid and pore compressibilities, k is the constant and uniform fault permeability, and η the fluid dynamic viscosity.
The fault is initially fully locked (zero slip rate) and the uniform initial fluid pressure field is equal to p o . We investigate sustained fluid injection for t > 0 under two different scenarios: either at constant volume rate or at constant overpressure. The solutions of both boundary value problems for the two-dimensional diffusion equation are well-known (Carslaw & Jaeger 1959) and can be written in the following functional form where ∆p * is a characteristic pressure and Π is the dimensionless injection-induced overpressure. These solutions of the diffusion equation (3) notably depend on the fault hydraulic diffusivity α = k/Sη [L 2 /T ] via the well-known diffusion characteristic length √ αt.
For a constant volume rate injection from a point source, the two-dimensional flow solution is given in polar coordinates by (section 10.4, eq. 5, Carslaw & Jaeger 1959): where Q w is the constant injection volume rate [L 3 /T ], w h is the fault thickness [L], and E 1 is the exponential integral function. Note that the product kw h is often denoted as the fault hydraulic transmissivity For the case of an injection from a finite wellbore at constant overpressure, the solution reads (section 13.5, eq. 6, Carslaw & Jaeger 1959): where ∆p w is the applied constant overpressure at the wellbore of radius r w , and J 0 and Y 0 are the zero-order Bessel functions of the first and second kind, respectively.
The uniform initial fluid pressure and stress fields must satisfy the condition |τ o | < f σ o on the fault plane, where σ o = σ o − p o is the initial effective normal stress. This condition means no activation of fault slip prior to the start of the injection. Fault slip starts when the fluid pressure increase is sufficient to reach the Mohr-Coulomb shear failure criterion. The ensuing aseismic rupture grows due to the direct effect of the fluid pressure that reduces locally the fault strength in Eq.
(2), and due to the quasi-static nonlocal elastic integral operator (1) that operates over the fault slip distribution δ and determines the local shear stress change consistent with the Mohr-Coulomb strength condition.

Numerical methods
We developed a fully implicit boundary-element-based solver with an elasto-plastic-like constitutive interfacial law to solve simultaneously the quasi-static elastic equilibrium (1) and the activation criterion (2) knowing the semi-analytical fluid pressure evolution on the fault (either Eq. (5) or (6) depending on the injection scenario). The main features of our numerical solver are briefly described below.

Spatial discretization
We discretize the fault plane Γ using an unstructured Delaunay triangulation Γ = ∪ N E k=1 Γ k , where Γ k is the k-th triangular element and N E is the total number of elements in the mesh (see Fig. 5a-b for examples of the unstructured meshes used). We use the displacement discontinuity method (Crouch & Starfield 1983) to solve the quasi-static elastic equilibrium and employ piecewise quadratic triangular boundary elements (Nikolskiy et al. 2013). The analytical integration of the traction kernel for a generic triangle with quadratic shape functions was obtained by Nikolskiy et al. (2013) for the three components of the displacement discontinuity. Our implementation considers boundary elements with six collocation points such that N C = 6N E is the total number of collocation points in the mesh, and N = 3N C is the total number of degrees of freedom.
In the configuration investigated here where the principal shear direction is known and does not change, the quasi-static elastic equilibrium, Eq. (1), can be written in the x-direction of the global reference system only. Nevertheless, we use a fully 3D collocation displacement discontinuity method and instead expressed it in the local reference system x k ; e k 1 , e k 2 , e k 3 of each k -th boundary element, where e k 1 and e k 2 are two unit vectors tangent to Γ k (and mutually orthogonal) and e k 3 is a unit vector normal to Γ k such that n + = −n − = −e 3 where "+" and "-" refer to the upper and bottom fault faces. The spatially discretized form of the quasi-static elastic equilibrium in the local reference system of the boundary elements is where t ∈ R N is the total traction vector, t o ∈ R N is the initial total traction vector, d ∈ R N is the displacement discontinuity vector, and E ∈ R N ×N is the collocation boundary element matrix which is dense and non-symmetric. We approximate this dense boundary element matrix E using a hierarchical matrix E H representation. Using a hierarchical matrix approximation allows to significantly reduce the memory requirements and speed up algebraic operations for an otherwise computationally expensive matrix (see Ciardo et al. 2020, and references therein for further details).
The system of equations (7) is arranged in order that t = t m i = (t 1 1 , t 1 2 , t 1 3 , t 2 1 ..., t N C 3 ), where the index i = 1, 2, 3 denotes the local components with regard to the local reference systems of the boundary elements, and the index m = 1, ..., N C denotes the collocation points.
Finally, the spatially discretized form of the Mohr-Coulomb shear failure criterion is solved locally (at the collocation point level) and is written as where τ m = (t m 1 , t m 2 ) is the local shear traction vector at the m-th collocation point, and t ,m 3 = t m 3 − p m is the normal component of the local effective traction vector at the m-th collocation point.

Time integration
We use a backward Euler time integration scheme. Let ∆X = X n+1 − X n be the increment of a generic variable X from the time t n to the time t n+1 with ∆t = t n+1 − t n being the time step. Taking the time derivative of Eq. (7) and using the definition of Terzaghi's effective stress such that t = t + p , where t ∈ R N is the effective traction vector and p ∈ R N is the fluid pressure vector, we arrive to the following fully discretized incremental form of the quasi-static elastic equilibrium Note that ∆p = (0, 0, ∆p 1 , 0, ..., ∆p N C ) as the fluid pressure affects only the normal components of the effective traction vector.
Eq. (9) is a nonlinear system of N equations for ∆d that is solved via a Newton-Raphson scheme. Note that ∆t depends nonlinearly on ∆d due to the Mohr-Coulomb shear failure criterion on the fault interface (8).

Integration of the constitutive interfacial law
The Mohr-Coulomb criterion (8) defines two different modes of contact at every collocation point which are either locked or in frictional sliding. On the contrary to the rigid-plastic approach developed in Ciardo et al. (2020), we solve here for ∆t (∆d) using an elasto-plastic constitutive interface relation between the local effective traction t and displacement discontinuity d. In other words, we regularize the fault frictional contact behavior by allowing a degree of elastic displacement discontinuity and write the increment of displacement discontinuity vector as the sum of an elastic part ∆d e and a plastic part ∆d p , ∆d = ∆d e + ∆d p .
We introduce a linear elastic relation between the increment of effective tractions ∆t and the elastic part of the displacement discontinuity, where D ∈ R 3×3 is a diagonal elastic stiffness matrix containing the shear and normal stiffness components of the local elastic springs (of dimension F/L 3 ) modeling the fault elastic response. Numerically, we consider sufficiently large values of the elastic stiffness components such that ∆d p ∆d e , whenever plastic flow occurs. Note that the minus sign in the above directly comes from our convention of signs for tractions (positive in compression) and displacement discontinuities (positive in opening).
The Mohr-Coulomb yield function in the local reference frame of a triangular displacement discontinuity element can be rewritten as If F < 0, the mode of contact is elastic and no frictional sliding occurs ∆d p = 0, whereas if F = 0, plastic frictional contact occurs and thus frictional sliding ∆d p = 0.
We use a non-associated Mohr-Coulomb flow rule with zero dilatancy to describe the evolution of the plastic part of the displacement discontinuity: where ∆γ ≥ 0 is the plastic multiplier increment and G is the plastic flow potential.
The previously mentioned inequalities for the yield function and plastic flow can be re-written as the Karush-Kuhn-Tucker conditions of elastoplasticity: For a given increment of the total displacement discontinuity vector ∆d, the system of Eqs. (10) to (14) has to be solved in order to obtain the plastic multiplier increment ∆γ, and consequently the increment of plastic fault slip ∆d p and the incremental change of effective tractions ∆t . We use a classical elastic predictor-plastic corrector algorithm to solve this system of equations. The elastic predictor-plastic corrector algorithm is a two-step procedure in which the two possible modes of contact, elastic and plastic, are solved sequentially and the final solution is chosen as the only one satisfiying the yield inequality (see for example de Souza Neto et al. 2008 for details). Note that this algorithm is executed locally at every collocation point for every Newton-Raphson iteration for a given increment of time. As usual in elastoplasticity, it is critical to use the consistent tangent operator during the Newton-Raphson iterations in order to achieve quadratic convergence. We recall the expression of this consistent tangent operator in the Appendix A.

Implementation details
Since we do not consider fluid flow being affected by mechanical deformation, the numerical time-stepping scheme consists only of solving a Newton-Raphson scheme for the mechanical equilibrium and the local elasto-plastic relation at time t n+1 knowing the fluid pressure increment given by the semi-analytical solutions (Eqs. (5) or (6)). Convergence of the non-linear Newton-Raphson solver is reached when the relative increment of the norm of the displacement discontinuity vector between two consecutive iterations falls below 10 −4 . The linear tangent mechanical system at each Newton step is solved using a biconjugate gradient stabilized iterative solver (BiCGSTAB) with a tolerance set to 10 −4 .
It is worth noting that the exponential integral function in Eq. (5) is log singular at the injection point r = 0. Hence, for the case of injection at constant volume rate, the fluid pressure is actually unbounded right at the origin. There is thus an evolving but rather small (comparing to the rupture size) circular region in the neighborhood of the injection point where the fluid pressure p is larger than the initial effective normal stress σ o . Even though our model does not account for fault opening, the effect of this singularity is negligible in our numerical solutions, since, even for relatively high resolution, p < σ o at every single discrete (collocation) point of the simulations.
3 Self-similar rupture growth due to injection at constant volume rate

Scaling and similarity
We first investigate the case of a constant volume rate injection, where the fluid pressure evolution driving the rupture is given by Eq. (5). Such a diffusion solution is self-similar: the pressure is only function of the self-similar variable r/ √ 4αt, where L(t) = √ 4αt is the characteristic diffusion lengthscale which corresponds to the evolution of the fluid pressure front disturbance. Because no other time or length scales enter the problem, the quasi-static rupture will also evolve in a self-similar fashion. This result is essential for the problem addressed in this section. We denote the a priori unknown rupture shape as R(t) = {(x, y) : δ(x, y, t) = 0} and scale it as where S is the dimensionless rupture front and R(t) the characteristic rupture lengthscale. Moreover, we define the amplification factor λ that relates the instantaneous rupture characteristic scale R(t) to the nominal location of the fluid pressure front L(t), such that R(t) = λL(t). A value of λ > 1 indicates that the rupture lengthscale is greater than the fluid pressure front radius, whereas a value of λ < 1 indicates the opposite.
We can scale the spatial variables (x, y) with the diffusion lengthscale L(t) (or alternatively with the characteristic rupture scale R(t)) while the characteristic fluid pressure scale is directly given by ∆p * in Eq. (5). Introducing these characteristic scales in the Mohr-Coulomb and elasticity equations, allows to close the scaling of the problem as: where the characteristic slip is given by δ c (t) = f ∆p * L(t)/µ (or alternatively by f ∆p * R(t)/µ if R(t) is used for the spatial scale). Using this scaling, the dimensionless form of the problem depends on only two dimensionless parameters: and the Poisson's ratio ν. The dimensionless rupture shape S thus depends on both ν and T . The parameter T is similar to the one found by Bhattacharya & Viesca (2019) in their 2D plane-strain model of a frictional shear crack driven by injection at constant pressure, whereas the second dimensionless parameter, the Poisson's ratio, arises from the three-dimensional nature of the problem considered here, in which the rupture propagates in mixed mode (II+III) with an a priori unknown shape.
The parameter T , named as fault stress parameter by Bhattacharya & Viesca (2019), is crucial in the present study. It encapsulates the information about the initial state of stress acting on the fault and the characteristic injection pressure. More specifically, the numerator of T , 1 − τ o /f σ o , is a measurement of the "distance" to failure under pre-injection ambient conditions (close to zero for a critically stressed fault, and close to one for a fault initially far away from frictional failure), whereas the denominator ∆p * /σ o is an overpressure ratio which measures the amount of pressurization due to injection with regard to the initial effective normal stress. Both numerator and denominator are indeed independent parameters of a more general model for a shear crack obeying slip-weakening friction in 2-D elastic media investigated by Garagash & Germanovich (2012).
The fault stress parameter varies between 0 and +∞. The limiting values of T are associated with end-member scenarios that are relatively similar to the ones identified by Garagash & Germanovich (2012) and Bhattacharya & Viesca (2019) for two-dimensional problems involving injections at constant pressure. For small values of T , the condition 1 − τ o /f σ o ∆p * /σ o must be satisfied. This means that the fault is "critically stressed" with regard to the overpressure ratio. For large values of T , the condition ∆p * /σ o 1 − τ o /f σ o must be satisfied, so that the fault is "marginally pressurized" with regard to the level of stress criticallity. Hence, following these prior studies, we denominate the corresponding end-member scenarios as a critically stressed fault (T 1) and a marginally pressurized fault (T 1).
It is worth noting that the characteristic pressure ∆p * = Q w η/4πkw h increases (and therefore the fault stress parameter T decreases) not only when the injection volume rate Q w grows, but also when the fluid viscosity η increases or the fault hydraulic transmissivity kw h decreases. On the other hand, large values of T might be eventually upper bounded if the fluid pressure near the injection point is high enough to make fault opening a significant mechanism driving the propagation of the rupture. Such problem has been already addressed in 2-D (Azad et al. 2017).

Analytical solution for circular ruptures
We first consider the particular case where the Poisson's ratio ν is set to zero such that the solution of the problem only depends on the value of T . In this limit (ν = 0), the rupture front R(t) is circular because the energy release rate for an axisymmetric shear load distributes uniformly along the circular crack front (see Appendix B). The dimensionless rupture shape S is thus the unit circle and R(t) is simply the rupture radius. For such a frictional shear crack with a constant friction coefficient, there is no fracture energy spent during propagation. The condition for quasi-static crack propagation then reads (see Appendix B) where the axisymmetric shear stress drop is given by The stress drop can be viewed as the contribution of two terms: a constant and negative term τ o − f σ o which is the difference between the initial shear stress τ o and the initial fault strength f σ o , and an axisymmetric and positive term f ∆p * Π(r, t) capturing the local reduction of fault strength due to fluid injection. After a change of variable s = r/R and incorporating the previous definition of the amplification factor λ = R(t)/L(t) into Eq. (18), the condition for quasi-static crack propagation can be rewritten in dimensionless form as As expected from scaling analysis, T is the only governing dimensionless parameter and, as a result, there is a unique value of λ for each value of T . The integral in Eq. (20) can be evaluated analytically to obtain the following implicit equation for λ as function of T where γ = 0.577216... is the Euler-Mascheroni's constant and 2 F 2 [ ] is the generalized hypergeometric function. The relation (21) is plotted in Fig. 2. This figure shows that a critically stressed fault (T 1) is characterized by a rupture that largely outpaces the fluid pressure front (λ 1). On the other hand, a marginally pressurized fault (T 1) is characterized by a rupture that substantially lags the fluid pressure front (λ 1).
The limiting behavior of Eq. (21) for small and large λ allows to obtain simple closed-form expressions for the corresponding end-member cases. In the limit of a critically stressed fault (λ 1), Eq. (21) can be asymptotically expanded as T ∼ 1/(2λ 2 ) + O(1/λ 4 ), leading to the following asymptotic solution for the amplification factor Likewise in the limit of a marginally pressurized fault (λ 1), Eq. (21) follows the asymptotic expansion Since in the critically stressed limit, the pressurized fault patch is small compared to the rupture area, the fluid pressure perturbation can be approximated by a monopole distribution. Such approximation is, in terms of local reduction of fault strength, equal to a point force given by f ∆p * Π(r, t) ≈ f ∆p * +∞ 0 Π(r, t)rdrδ dirac (r)/r = 2αf ∆p * tδ dirac (r)/r, where δ dirac is the Dirac delta function in polar coordinates. Replacing this approximation in Eq. (19), and then evaluating the crack propagation condition, Eq. (18), leads equivalently to the asymptotic solution for λ in the critically stressed limit, Eq. (22).
On the other hand, in the marginally pressurized limit where the crack size is small compared to the pressurized area, the fluid pressure perturbation within the crack can be approximated by considering the behavior of the the exponential integral function in Eq. (5) for small values of its argument. Such approximate injection-induced local reduction of fault strength is f ∆p * Π(r, t) ≈ −f ∆p * 2 ln r/ √ 4αt + γ . Again, replacing this approximation in the crack propagation condition, Eq. (18), leads alternatively to the asymptotic solution for λ in the marginally pressurized limit, Eq. (23).
The purpose of the previous analysis is to highlight the asymptotic form of the driving forces related to the two end-member cases. As we will see later, both fault slip and rupture shape will also show well-defined asymptotic behaviors, and thus the results can be directly associated with the type of force that drives the crack growth in both limiting cases.
The asymptotic solutions (22) and (23) are shown in Fig. 2 together with the general solution given by Eq. (21). Note that the transition between both propagation regimes (defined as λ = 1) occurs at T ≈ 0.5915.
Another interesting analytical result is the rupture speed V r that decreases with the squared root of time and is simply given by (marginally pressurized fault) (critically stressed fault) The singularity of the rupture speed at t = 0 is a consequence of the self-similar diffusion lengthscale √ 4αt and the abscense of inertia. In addition, the (de)acceleration of the rupture front is equal to −λ √ α/(2t 3/2 ). This power-law (de)acceleration is comparable to tensile hydraulic fracture propagation under a constant injection rate albeit at a different power-law of time (Detournay 2016).

Numerical solution for circular ruptures
The numerical solution for circular ruptures allows us to obtain the axisymmetric self-similar slip profiles and also to verify our numerical solver against the previously derived solution for the amplification factor λ(T ). We use the boundary-element-based solver previously described in section 2.2 and run seven simulations for values of T = 0.001, 0.01, 0.1, 0.7, 2.0, 4.0 and 7.0 with ν = 0. We perform 10 fully implicit time steps for each simulation. Owing to the self-similarity of the problem, the slip profiles at different times collapsed into one single curve under the scaling of Eq. (16). As a consequence, there is a unique dimensionless slip profile for a given value of the fault stress parameter T . The unique self-similar slip profiles for the different values of T are shown in Figs. 3b and 3d for critically stressed and marginally pressurized cases, respectively. Moreover, in Appendix C, we derive closed-form analytical expressions for the self-similar slip profiles in the limiting cases of critically stressed (λ 1) and marginally pressurized (λ 1) faults, that are:

Axisymmetric slip profiles and accumulated fault slip at the injection point
in the marginally pressurized limit, and in the critically stressed limit, wherer = r/R(t) is the self-similar radial coordinate. Eq. (26) is indeed valid for r L(t) and it corresponds to the "outer" solution in the critically stressed limit. Both analytical expressions are shown in Figs. 3b and 3d together with the numerical results.
The accumulated fault slip at the injection point, δ(r = 0, t), is plotted in Fig. 7c as a function of the fault stress parameter T . Note that δ(r = 0, t) is normalized by the position of the fluid pressure front L(t) on the left axis and by the rupture radius R(t) on the right axis. δ(r = 0, t) is easily obtained from Eq. (25) in the marginally pressurized limit as whereas in the critically stressed limit, The prefactor in the previous equation is obtained numerically and it is approximately 3.5 (see Fig. 7c).
Eqs. (25) to (28) confirm that the relevant scale for the shear stress is f ∆p * , as chosen in the scaling analysis. Also, it becomes now clear that the relevant lengthscale in the problem depends on the limiting case under consideration; for marginally pressurized faults, the relevant lengthscale is the rupture radius R(t), whereas for critically stressed faults, the proper lengthscale is the nominal radius of the pressurized fault patch L(t).

Rupture radius and solver verification
In order to determine numerically the instantaneous rupture radius R(t) at every time step t n and verify our numerical solver against the analytical solution, Eq. (21), we solve numerically for R(t n ) by searching for the position of zero slip δ(R, t n ) = 0. As the solution for slip is axisymmetric in the limit of ν = 0, we search in fact for the zeros along the entire rupture front by taking 100 equally-spaced values of the angular cylindrical coordinate θ ∈ [0, 2π). In this way, we build the rupture front and compute finally the instantaneous rupture radius by fitting the equation of a circle centered at the origin to the zeros found for all the values of θ considered. Fig. 4a shows the results for the rupture radius as a function of the nominal location of the fluid pressure front L(t) = √ 4αt for different values of the fault stress parameter T . In such plot, self-similar solutions for the rupture growth in the form R(t) = λL(t) are represented by straight lines that cross the origin and have a slope equal to the amplification factor λ. We estimate numerically the amplification factor λ for each value of the fault stress parameter T , by simply averaging the ratios R(t n )/L(t n ) over the different time steps of the simulation.
The numerical results for the amplification factor λ are displayed in Fig. 4b together with the analytical solution, Eq. (21). The numerical results are in excellent agreement with the theoretical predictions. This plus the previous comparison with the asymptotics of fault slip allows us to verify our numerical solver before exploring the case of non-circular ruptures, which is solved by numerical means only. The relative error between the numerical results and the exact analytical solution for the amplification factor λ is showed in the inset of Fig. 4b and is approximately below 1%.

Numerical solution for non-circular ruptures
We move now to the more general case where the Poisson's ratio is different than zero. It is important to recall that the rupture shape R(t) is not known a priori but, of course, remains self-similar. In order to cover the relevant parameter space, we run 21 simulations for the same seven values of the fault stress parameter T considered in the previous section, 0.001, 0.01, 0.1, 0.7, 2.0, 4.0 and 7.0, for three values of the Poisson's ratio ν = 0.15, 0.30, and 0.45.

Rupture shape
We quickly recognize in our simulations that the ruptures evolve systematically in a nearly elliptical shape. We also observe that the aspect ratio of the ruptures depends strongly not only on the Poisson's ratio but also on the fault stress parameter T , i.e., on the initial stress state and the driving force itself. Snapshots of two rupture simulations having the same Poisson's ratio but different values of the fault stress parameter are shown in Figs. 5a-b. It is clear that the aspect ratio for critically stressed faults (Fig. 5a, T = 0.001, λ 1) is higher than the aspect ratio for marginally pressurized faults (Fig. 5b, In order to quantify the rupture shape, we perform a nonlinear regression of the rupture front at every time step assuming an elliptical shape: where a and b are the semi-major and semi-minor axes of the ellipse. We use the same procedure described previously to estimate the rupture front, with the only difference that now the spatiotemporal evolution of slip is no longer axisymmetric. Typical ellipsoidal fits of the rupture front are displayed in Figs. 5a-b. The results for the aspect ratio a/b as a function of the fault stress parameter T and the Poisson's ratio ν are summarized in Fig. 5c. Note that the aspect ratio is time-invariant due to the self-similarity of the problem, thus, we average the aspect ratio over the simulations' time steps for better accuracy.  (d) Aspect ratio a/b with increased resolution for the Poisson's ratio ν for the two end-member cases, T = 0.001 (λ ≈ 22.37, critically stressed faults) and T = 7.0 (λ ≈ 0.03, marginally pressurized faults).
displays clearly two asymptotic behaviors of the aspect ratio for the two end-member cases of a critically stressed fault and a marginally pressurized fault. Additional simulations were run to better explore the dependence on Poisson's ratio for T = 0.001 (λ ≈ 22.37, critically stressed faults) and T = 7.0 (λ ≈ 0.03, marginally pressurized faults). The results are shown in Fig. 5d. We notably found by numerical observation that the aspect ratio grows asymptotically with the Poisson's ratio as It is interesting to note that the aspect ratio for critically stressed faults, a/b ∼ 1/(1 − ν), is similar to the results obtained by Gao (1988) and Elie et al. (2006) for three dimensional planar shear-crack under uniform remote loading and a uniform energy release rate along the crack front.

Generalized amplification factor λ
Since the rupture front R(t) is self-similar and nearly elliptical, we can write a(t) = λ a L(t) and b(t) = λ b L(t), where λ a and λ b are the corresponding amplification factors for the semi-major and semi-minor axes of the rupture front, respectively. Note that λ a and λ b depend only on the fault stress parameter T and the Poisson's ratio ν. The evolution of a(t) and b(t) as a function of the fluid pressure front location L(t) = √ 4αt is shown in Fig. 6a for different values of the fault stress parameter T and the Poisson's ratio ν. In this figure, we also include the reference circular case ν = 0, and indicate the meaning of λ a , λ b , and λ as the slopes of the straight lines for a(t), b(t), and R(t) (of the reference circular case), respectively. Fig. 6a shows that the major and minor axes for the elliptic ruptures (ν = 0) lie about the radius for the reference circular solution (ν = 0). For a given value of T and any value of ν, we find that the geometric mean of a(t) and b(t), a(t)b(t), is equal to the radius R(t) of the circular crack solution for the same value of T (and ν = 0). This equality is equivalent to the equality of amplification factors This is demonstrated in the inset of Fig. 6a in which the numerical values of λ num = √ λ a λ b for all values of T and ν are plotted against the exact solution for circular ruptures λ circ .
The numerical results for Eq. (31) are plotted in Fig. 6b together with the analytical solution, Eq. (21), for all values of T and ν. In the inset, the relative difference between the numerical results and the analytical solution are also displayed. Fig. 6b thus demonstrates that Eq. (31) is a generalization of the amplification factor that is now valid for any value of the Poisson's ratio. In the particular case of ν = 0, Eq. (31) reduces simply to λ = R(t)/L(t), as originally defined when deriving the analytical solution for the circular rupture case.

Poisson's ratio-independent rupture area
The generalized amplification factor λ = √ λ a λ b has a clear physical meaning. Indeed, it is equivalent to the squared root of the ratio between the instantaneous elliptic rupture area A r (t) = πa(t)b(t) and the instantaneous pressurized area πL 2 (t): λ = A r (t)/(πL 2 (t)) such that the evolution of the rupture area A r (t) is simply given by and is thus independent of the value of the Poisson's ratio ν. The Poisson's ratio (together with the value of T ) modifies the shape of the ruptures, which are more or less elongated, but it does not modify the rupture area, which solely depends on T . The rupture area A r evolves linearly with time and proportionally to the injected volume (∝ V ) for such a constant injection rate case.
Furthermore, for the two end-member cases of critically stressed and marginally pressurized faults, Eqs. (22) and (23) lead to simple closed-form expressions for the evolution of the rupture area as function of the fault stress parameter T A r (t) ∼ 2παt/T for critically stressed faults (T 1, λ 1) A r (t) ∼ παe 2−γ−T t for marginally pressurized faults (T 1, λ 1) It is worth noting that David & Lazarus (2021) obtained recently a somewhat similar result in their study of tensile crack growth under the Paris' fatigue law (with a uniform energy release rate being a limiting case). They found, also by numerical observation, that a circular crack solution is sufficient to predict the area of a rupture for any non-circular crack.

Rupture front R(t) for the end-member cases
The asymptotic expressions for the aspect ratio (30) plus Eq. (33), lead to the following closed-form expressions for the evolution of the semi-major a(t) and semi-minor b(t) axes of the rupture front for critically stressed faults (λ 1): and for marginally pressurized faults (λ 1) Since the rupture front is quasi-elliptical, Eqs. (34) and (35) fully define the spatiotemporal evolution of the rupture front R(t) for the end-member cases.

Non-axisymmetric slip profiles and accumulated fault slip at the injection point
The non-axisymmetric self-similar slip profiles are unique for a given combination of T and ν. Some typical slip profiles along the x -axis (normalized by a(t)) are shown in Figs. 7a-b for different values of the fault stress parameter T and ν = 0.3. The accumulated fault slip at the injection point δ(r = 0, t) is plotted in Fig. 7c for all simulations as a function of the fault stress parameter T and the Poisson's ratio ν. In Fig. 7c, we include the circular rupture case (ν = 0) and δ(r = 0, t) is further normalized in Fig. 7d by the geometric mean a(t)b(t) which is Poisson's ratio-independent and in the limit of ν → 0 corresponds to the rupture radius R(t).
Figs. 7c-d shows that the accumulated fault slip at the injection point decreases for increasing values of the Poisson's ratio. In addition, we recover a similar scaling for δ(r = 0, t) that in the circular rupture case: δ(r = 0, t) ∼ f ∆p * R(t)/µ in the marginally pressurized limit, with the characteristic rupture scale R(t) taken as a(t)b(t) in Fig. 7d; and δ(r = 0, t) ∼ f ∆p * L(t)/µ in the critically stressed limit.

Aseismic moment
The scalar moment release M 0 is a key measurement in seismology to quantify the potency of a rupture (Aki & Richards 2002). In our planar fault model with a uniform direction of slip δ, the time-dependent aseismic moment (we focus on the case of a circular rupture) is M 0 (t) = 2πµ δ(r, t)rdr, where R(t) is the evolving rupture radius. We can thus compute the aseismic moment numerically for all the seven values of T considered. Furthermore, we can use the asymptotic solutions of the slip distribution, Eqs. (25) and (26), to derive closed-form expressions for the limiting behaviors of the aseismic moment. We obtain that in the marginally pressurized limit, and in the critically stressed limit.
Both previous equations provide the proper scaling of the aseismic moment. We use these scalings to normalized the numerical results that are presented in Fig. 8 as a function of the fault stress parameter T . In this figure, we include the corresponding prefactors of Eqs. (36) and (37). Note that the prefactor 8/3 in the critically stressed limit is in good agreement with the numerical solution, despite the slip profile being approximated by the "outer" asymptotic solution of this limit only.
4 Non-self-similar rupture growth due to injection at constant pressure

Scaling and similarity
Under the scenario of constant pressure injection, the evolution of fluid pressure given by Eq. (6) is no longer self-similar. This is due to the presence of a finite wellbore radius r w where the pressure is set constant which introduces a new characteristic length in the problem. As a result, the frictional rupture will not evolve self-similarly, like for a constant injection rate. We recall that self-similar solutions always correspond to idealized models in which the dimensional parameters of the independent variables (space and time in our case) are equal to zero or infinity (Barenblatt 1996). The infinitesimal nature of the constant-volume-rate fluid source of the previous section and its subsequent self-similarity, can be seen in fact as an intermediate-asymptotic behavior of a more realistic physical system in which both the fluid source and the domain are finite. In this view, the solution of section 3 is valid for times t r 2 * /α, where r * is the characteristic length of the fluid source, and for t R 2 * /α, where R * is the characteristic length of a finite domain. Note that the introduction of, for instance, a frictional lengthscale in the self-similar problem of injection at constant volume rate, would also cause the loss of self-similarity in the model.
The scaling thus differs slightly from the scaling of the previous section. The finite wellbore radius r w introduces a characteristic diffusion timescale t c = r 2 w /α with α the fault hydraulic diffusivity, such that we obtain: where ∆p w is the constant overpressure imposed at the wellbore, and δ c = f ∆p w r w /µ is the characteristic slip.
As already mentioned, the loss of self-similarity is due to the finite size of the wellbore. In fact, radial flow from an infinitesimal fluid source at constant pressure is not physically possible. Injection of a finite volume from an infinitesimal fluid source in such geometrical conditions always leads to infinite pressure.
The dimensionless solution of the problem depends now on three dimensionless parameters, a slightly different fault stress parameter T that is a function of the constant wellbore overpressure ∆p w the Poisson's ratio ν, and the dimensionless time αt/r 2 w . The limiting values of T are determined by the condition for fault slip activation, f ∆p w ≥ f σ o − τ o , and the condition for no slip prior to injection, f σ o − τ o > 0. Together, these conditions imply that T varies between 0 and 1, so that the fault stress parameter is now upper bounded, unlike the case of injection at constant volume rate in which T can theoretically go up to +∞.
The limit of T → 1 is characterized by the condition f ∆p w → f σ o −τ o . This condition can be interpreted as a scenario in which the pressure at the fluid source, ∆p w , is just enough to activate fault slip. This is the reason why such end-member case has been named in prior studies as marginally pressurized faults (Bhattacharya & Viesca 2019, Garagash & Germanovich 2012. On the other hand, considering that ∆p w is always positive and finite, the limit of T → 0 is associated with the condition τ 0 → f σ 0 . This condition represents the case of faults that are about to fail before injection, and is thus named as critically stressed faults. Unlike the problem of injection at constant volume rate in which the fluid pressure near the injection point is always increasing, here the pressure at the wellbore is fixed and thus the terminology of critically stressed and marginally pressurized faults is unambiguous. Note that σ 0 > ∆p w > 0, with the upper bound being the transition to fault opening that we do not cross in this study.

Approximate analytical solution for circular ruptures
The solution of the diffusion equation for injection at constant pressure from a finite source, Eq. (6), does not allow to treat the problem of a circular rupture analytically. Nevertheless, as we are interested in solutions for times that are in general large compared to the characteristic diffusion time, t c = r 2 w /α, we can use the following asymptotic expansion for the instantaneous pressure profile that is valid for t t c (Jaeger 1955) wherer = r/r w andt = αt/r 2 w are the dimensionless radial coordinate and the dimensionless time, respectively, and γ = 0.577216... is the Euler-Mascheroni's constant. Fig. 9a shows the comparison between the exact numerical solution for the fluid pressure profile, Eq. (6) (solved via numerical inversion of the Laplace transform, Stehfest 1970) and the asymptotic expansion (41). In Fig. 9a, we also include an asymptotic expansion of Π(r, t) with higher order corrective terms that are function ofr and the successive time derivatives of the terms in curly brackets in (41) (Jaeger 1955). The higher order terms have cumbersome expressions but are necessary to capture the "near front" behavior of the fluid pressure profile as shown in Fig. 9a. However, for the sake of simplicity, we neglect these corrective terms in the following.
Similarly to the case of injection at constant volume rate, we define the instantaneous rupture radius R(t) and use the condition for quasi-static propagation of a circular crack with zero fracture energy under axisymmetric shear load, Eq. (18), with now ∆τ (r, t) = τ o − f (σ o − ∆p w Π(r, t)). After some algebraic operations, this propagation condition can be rewritten as where T is the fault stress parameter defined in Eq. (40) for the constant pressure injection case. We approximate Π(r, t) = (p(r, t) − p 0 ) /∆p w by the asymptotic expansion (41). Note that in Eq. (41), one could consider to drop the term of O 1/ ln(t) 2 and use a first order approximation for Π(r, t) instead; however, we found that better results are systematically obtained by keeping the second order term in any further mathematical operation and performing first order approximations afterwards.
The integration limits of Eq. (42) have to be considered carefully, since the asymptotic expansion for Π(r, t) gives non-physical values that are greater than unity for r/r w < 1 and negative for r/r w beyond the intersection with the abscissa (see Fig. 9a). In fact, the intersection with the abscissa defines conveniently a nominal position of the fluid pressure front,L(t)/r w = exp(−(1/2) (2γ − ln (4t)) 2 / (3γ − ln (4t))), that is given at the first order byL where c 1 = e ln(4)−γ = 2.245838... ≈ 2.25. The position of the fluid pressure frontL(t) given by Eq. (43) is shown in Fig. 9a at different dimensionless times. With a change of variable s = r/R and taking care of the integration limits as discussed before, we can rewrite Eq. (42) in dimensionless form as where β 0 = r w /R, and β = 1 if R ≤L, or β =L/R otherwise (R >L).
Eq. (44) can be solved to obtain the evolution of the normalized rupture radius R(t)/r w as a function of the dimensionless time αt/r 2 w and the fault stress parameter T . The solution is piecewise due to the piecewise definition of β that indeed separates the two possible rupture regimes. One regime is characterized by R(t) <L(t) which represents a rupture front that lags the fluid pressure front, whereas the other regime is characterized by R(t) >L(t) in which the rupture front outpaces the fluid pressure front.
In addition, the first integral of the left-hand side of Eq. (44) can be neglected if we assume that the rupture radius R(t) is much bigger than the wellbore radius r w , so that β 0 = r w /R 1. Hereafter, we consider β 0 = 0. Let us first consider the case in which R(t) <L(t). After integrating Eq. (44) with β = 1, we obtain at the first order the following explicit expression for the evolution of the normalized rupture radius in the form of a power-law where g (T ) = c 2 e −c3T , with c 2 = e (2−γ)/2 = 2.036825... and c 3 = (ln (4) − γ) /2 = 0.404539....
Note that the transition between the two rupture regimes happens at a certain time t * when R(t * ) = L(t * ). Using the first-order Eqs. (43) and (45), we obtain that this transition time is

Finally, the solution for the case R(t) >L(t) is obtained by integrating Eq. (44) with β(t) =L(t)/R(t).
The solution for the rupture radius is now implicit and it is given by where f 1 (β) = 1 − 1 − β 2 , and f 2 (β) = ln 2β Eqs. (45) to (48) can be used to define a time-dependent amplification factor in the form λ(t) = R(t)/L(t).
Such approximate analytical solution for λ(t) is plotted in Fig. 9b at different dimensionless times, as a function of the fault stress parameter T .

Numerical solution for circular ruptures
We now solve numerically for the case of circular ruptures to obtain the evolution of the axisymmetric slip profiles δ(r, t). In addition, the computation of the slip profiles allows us to calculate numerically the rupture radius R(t) and test the accuracy of the approximate analytical solution derived in the previous section. For this purpose, we run six simulations for values of the fault stress parameter T = 0.01, 0.1, 0.3, 0.5, 0.7, and 0.9 and set ν = 0. We perform 9 fully implicit time steps per simulation for values of the dimensionless time logarithmically spaced between 1 to 10 8 . (a) Figure 9: (a) Instantaneous spatial profile of fluid pressure for injection at constant pressure, Π(r, t), as function of the dimensionless radial coordinate r/r w at different dimensionless times αt/r 2 w . Comparison between the exact solution given by Eq. (6), the asymptotic expansion (41) valid for αt/r 2 w 1, and an asymptotic expansion with higher order corrective terms (Jaeger 1955). (b) Approximate analytical solution for circular ruptures driven by injection at constant pressure, given by the amplification factor λ(t) = R(t)/L(t) as a function of the fault stress parameter T at different dimensionless times αt/r 2 w .

Rupture radius and comparison with approximate analytical solution
Based on the numerical solution of the axisymmetric slip profiles δ(r, t), we compute the instantaneous rupture radius R(t) at every time step for each simulation. We use the same procedure described in section 3.3 for building the rupture front and computing the rupture radius. The results are plotted in Fig. 11a together with the approximate analytical solution derived in the previous section.
The approximate analytical solution (valid for large times, αt/r 2 w 1) is in good agreement with the numerical results for values of T ranging from 0.1 to 0.7, with an average relative difference of about 5%. Near the limit of a marginally pressurized fault (T = 0.9), the analytical solution is less accurate (average relative difference around 20%) due to the fact that the assumption R(t) r w is not properly satisfied. On the other hand, near the limit of a critically stressed fault (T = 0.01), the analytical solution loses accuracy possibly due to the fact that the "near front" behavior of the fluid pressure profile is not well captured by the asymptotic expansion, Eq. (41). The absolute value of the relative difference between the approximate analytical solution and the numerical results is around 30% in average. Fig. 11b displays the numerical results for the time-dependent amplification factor λ(t) = R(t)/L(t) and the corresponding approximate analytical solution for it. Note that for values of T 0.7, the rupture lags the fluid pressure front at all times, whereas for T 0.01 the rupture outpaces the fluid pressure front at all times considered here. The case of intermediate values, 0.5 T 0.1, is interesting because a transition of propagation regime occurs at early times. Similarly to the case of injection at constant volume rate, the simulations show that ruptures evolve systematically in a nearly elliptical shape. We thus define the rupture front R(t) according to the equation of an ellipse (Eq. (29)), and compute the semi-major a(t) and semi-minor b(t) axes of the elliptical front following the same procedure described in section 3.4. Fig. 12a shows the temporal    12b displays the ratio between the geometric mean λ a (t)λ b (t) for ν = 0.3, where λ a and λ b are defined as λ a (t) = a(t)/L(t) and λ b (t) = b(t)/L(t), respectively, and the numerical values of the amplification factor λ(t) for the circular rupture case (ν = 0). We observe that like the case of injection at constant volume rate, the (now time-dependent) geometric mean √ λ a λ b is almost equal to the amplification factor λ(t) for circular ruptures, meaning that the rupture areas for the values of ν = 0 and ν = 0.3 are approximately the same for the values of T considered here.

Comparison between 2-D and 3-D models
We examine here the differences between our 3-D model and its counterpart in 2-D. In the twodimensional case, the diffusion of fluid pressure along the one-dimensional frictional interface is selfsimilar under both injection scenarios (constant volume rate and constant pressure) when considering a fluid point source (Carslaw & Jaeger 1959). Injection-induced fault slip will thus evolve in a self-similar fashion in both cases owing to the abscense of other lengthscales in the model.
The solution in 2-D elasticity for the evolution of the crack length under injection at constant pressure was presented by Bhattacharya & Viesca (2019). They showed that the amplification factor λ = (t)/ d (t), where (t) is the position of the crack tip (equal to the half-crack length) and d (t) is a nominal position of the fluid pressure front, is time-invariant and depends only on the fault stress parameter T . Interestingly, we found qualitatively the same response in our 3-D model but for injection at constant volume rate.
On the other hand, the solution of the 2-D model for injection at constant volume rate has not been presented yet. We derive such solution in the Appendix D and found that the amplification factor λ is time-dependent and follows exp −λ 2 /2 1 + λ 2 I 0 λ 2 /2 + λ 2 I 1 λ 2 /2 − 2λ/ where d = √ 4αt and x c = (f σ o − τ o ) / (f q w η/k). I 0 and I 1 are the Bessel functions of the first kind of zero and first order, respectively. Eq. (49) represents a unique relation between the amplification factor λ and the ratio between the position of the fluid pressure front d and the characteristic length x c , which is plotted in Fig. 13 together with the corresponding asymptotes for small and large λ (details in Appendix D).
The characteristic length x c corresponds to the position of the fluid pressure front at the onset of crack growth (activation of slip). Upon crack initiation, the rupture lags the fluid pressure front and expands faster than the diffusion of fluid pressure. The crack catches the fluid pressure front (λ = 1) when the normalized fluid pressure front¯ d ≈ 3.1435 or, in other words, when the fluid pressure front d has grown approximately three times the size x c necessary for the crack to start growing. After that, the rupture outpaces the fluid pressure front and the crack keeps propagating faster than the diffusion of fluid pressure until it reaches a steady propagation regime that is characterized by a constant rupture speed V R equal to (see Appendix D) where q w [L/T ] is the constant injection volume rate per unit fault thickness w h and unit out-of-the-plane length b, such that q w = Q w /w h b with Q w the injection volume rate L 3 /T . This response of the 2-D model under constant rate of injection has no analog in 3-D. Injection at constant volume rate in the 3-D model leads to a rupture speed that decreases with the squared root of time, (24)). Moreover, the relative position of the rupture front and the fluid pressure front is time-invariant in 3-D.
Our analysis shows that the response of the 2-D and 3-D models under the same injection scenario are different even qualitatively. These differences have to be carefully considered when linking theoretical and numerical predictions to laboratory measurements and field observations in which, generally, 3-D models prevail.

Assumption of constant friction
In the context of rock friction and earthquake mechanics, laboratory-derived friction laws (Dieterich 1979, Ruina 1983, Marone 1998) have been widely used to describe the entire spectrum of slip rates in natural faults (Scholz 2019). These empirical friction laws capture the dependence of friction on slip rate and the history of slidding (via a state variable) as observed during velocity-step laboratory experiments on bare rock surfaces and simulated fault gouge (Marone 1998). It seems then interesting to discuss under what conditions a constant friction coefficient can mimic more complex models of fluid-driven frictional ruptures with rate-and-state friction.
Results in 2-D antiplane elasticity of fluid-driven fault slip propagating on a strengthening (aging) rateand-state frictional interface have been notably reported by Dublanchet (2019), for the case of injection of fluid at a constant volume rate. Two distinct regimes of crack propagation were observed in Dublanchet's numerical simulations (following a first initial phase of slip rate acceleration): a stable crack growth that tends ultimately to a constant rupture speed, and an unstable crack growth in which the rupture speed blows up in a finite time. He observed that what determines the stable/unstable fault response is the sign of the difference between the residual shear stress left within the crack τ r and the initial shear stress resolved on the fault plane τ o . If τ r − τ o > 0 (a condition that is guaranteed in a constant-friction model), crack propagation is always stable, whereas if τ r − τ o < 0, the crack may evolve towards an instability. Such ultimate stability condition can be indeed understood under the classic Griffith energy-balance and the small scale yielding approximation (Rice 1968), as performed by Garagash & Germanovich (2012) for a fluid-driven slip-weakening frictional shear crack in 2-D.
The stable regime of crack propagation found by Dublanchet (2019) is indeed the most relevant in the context of aseismic ruptures. Furthermore, during such stable propagation regime, Dublanchet noticed that the crack behaves exactly as if it were governed by a constant friction coefficient within the slipping patch. This is because in that regime, the leading-order terms of the quasi-static elastic equilibrium are the nonlocal stress transfer along the fault and the effect of fluid pressure change on reducing the constant part of the shear strength in the rate-and-state friction law. It is not surprising then that our 2-D model of a constant-friction shear crack derived in the Appendix D and summarized in Fig. 13, shows qualitatively the same response (under the same injection scenario), notably the ultimate steady crack propagation regime characterized by a constant rupture speed V R . In our 2-D model, (50)), which depends on the constant friction coefficient f , hydraulic diffusivity α, injection rate q w , and residual τ r and initial τ o shear stress, in the same form as in the rate-and-state model (Eq. 23 in Dublanchet 2019), considering that f σ o is the residual shear strength in the constant-friction model.
A similar result, also in 2-D elasticity, has been recently obtained by Garagash (2021) using a different approach. He developed a Griffith-energy-balance-like equation of motion for the evolution of crack length on rate-and-state faults that he then applied to the study of slip transients due to point-force-like fluid injections. He notably showed that for injection at constant volume rate on neutraly and understress (with regard to the ambient slip rate) strengthening rate-and-state faults (obeying the slip law), the frictional ruptures expand initially within the limits of the pressurized fault patch and move faster than the latter, until they eventually outpaces the fluid pressure front and reaches also a terminal steady propagation regime characterized by a constant rupture speed. This response is again qualitatively the same of our 2-D model of a constant-friction shear crack, and the same found by Dublanchet (2019) for the stable propagation regime. Moreover, as pointed out by Garagash (2021), the ultimate behavior of the crack under these conditions is due to the diminishing effect of rate-and-state fracture energy in the Griffith energy balance compared to the effect of the fluid injection (in the energy release rate) with increasing rupture size, such that, in the limit of large-run-out rupture, the crack behaves as having zero toughness or, in other words, a friction coefficient that is constant.
Our discussion in 2-D elasticity suggests that the assumption of a constant friction coefficient describes to first order the behavior of rate-and-state friction under conditions that are relevant for the study of fluid-driven aseismic ruptures (rate-strengthening in both aging and slip laws, and approaching the large-run-out rupture regime), in which the frictional fracture energy can be neglected. We expect our results in 3-D to provide also first-order descriptions of fluid-driven aseismic ruptures in the context of rate-and-state friction, yet this assumption remains to be confirmed in future studies.

Implications to injection-induced seismicity
As suggested by a number of experimental and observational studies (Hamilton & Meehan 1971, Scotti & Cornet 1994, Bourouis & Bernard 2007, Guglielmi et al. 2015, Wei et al. 2015, Chen et al. 2017, Duboeuf et al. 2017, Eyre et al. 2019, Cappa et al. 2019, aseismic slip seems to be a frequent result of fluid injections into the subsurface and might play a significant role in injection-induced seismicity related to hydrocarbon and geothermal operations. It is thought that fluid motion drives aseismic slip which in turn transmits solid stresses that trigger part of the observed induced seismicity (Guglielmi et al. 2015, Wei et al. 2015, Duboeuf et al. 2017, Eyre et al. 2019, Bhattacharya & Viesca 2019. Our results open the possibility of quantifying this triggering mechanism in a three-dimensional scenario that is more realistic than previous two-dimensional models for fluid-driven aseismic fault slip (Eyre et al. 2019, Bhattacharya & Viesca 2019, Dublanchet 2019, Garagash 2021. First, we note that injection protocols generally consist of a series of injections conducted at a constant volume rate. Results of section 3 are thus more relevant for geo-energy applications. In particular, we showed that aseismic slip induced by injection at constant volume rate is self-similar in a diffusive manner. As a consequence, the rupture front expands proportionaly to the square root of time in the same way as the fluid pressure front does. Induced seismicity that is commonly considered to be driven by the direct effect of fluid pressure increase due to the square-root-of-time feature of seismicity clouds (Shapiro et al. 1997, might instead be controlled by the stress transfer of a propagating aseismic rupture (Bhattacharya & Viesca 2019). This would notably be the case of critically stressed fractures/faults in which the rupture front is predicted to be systematically ahead of the fluid pressure front. Our results suggest that assessing whether seismicity is induced by aseismic-slip stress transfer or fluid pressure increase might not be possible from the observation of square root time dependence of seismicity front alone.
Our model is, of course, idealistic in the sense that it represents a single and isolated fracture/fault in 3-D. Nonetheless, recent two-dimensional simulations of fluid induced aseismic slip in fractured rock masses have shown that the same patterns predicted by a single fracture in 2-D emerge collectively for a set of fractures (Ciardo & Lecampion 2021). Notably, a collective aseismic slip front outpaces the migration of fluids when the fracture network is in the critically stressed regime in a global sense (Ciardo & Lecampion 2021). In addition, field observations indicate critically stressed fractures/faults are likely to be preferred, high-permeability, conduits of fluid flow than the fractures/faults that are not optimally oriented with regard to the stress field (Barton et al. 1995). Together, these observations suggest that seismicity triggered by injection-induced aseismic slip might be indeed a general feature of reservoir rocks in response to fluid injections.
Moreover, if aseismic slip is the dominant mechanism for the triggering of seismicity, current estimates of reservoir hydraulic diffusivity α based on the spatio-temporal seismicity patterns (Shapiro et al. 1997) might be rather related to the quantity αλ 2 (see Eq. (32)), with λ being an equivalent amplification factor of the fractured rock mass. Such amplification factor would be intrinsically dependent not only on hydraulic properties of the fracture network, but also on the distribution of fracture orientations with regard to the stress field and the rate of injection.
Another finding of our study is related to the scaling relation between the aseismic moment release M 0 and the accumulated injected volume of fluid V . This type of relation has been extensively sought with the purpose of constraining the magnitude of injection-induced earthquakes based on operational parameters (McGarr 2014, van der Elst et al. 2016, Galis et al. 2017, McGarr & Barbour 2018. We found that the aseismic moment scales to the injected volume of fluid as M 0 ∝ V 3/2 , for injection at constant volume rate. Interestingly, the same power-law scaling has been found for self-arrested injection-induced seismic ruptures based on fully-dynamic rupture simulations and fracture mechanics arguments (Galis et al. 2017). We emphasize that our scaling relation is derived for purely aseismic (quasi-static) ruptures, whereas the relation found by Galis et al. (2017) represents seismic (dynamic) events.

Implications to seismic swarms, aftershock sequences and slow slip events
Seismic swarms and aftershock sequences are often characterized by diffusive spatiotemporal patterns that are thought to be caused by naturally injected fluids into fault zones (Bosl & Nur 2002, Miller et al. 2004, Parotidis et al. 2005, Chen et al. 2012, Hainzl et al. 2016, Ross et al. 2017, 2020. Moreover, natural fluid releases are likely represented by the sudden increase of injection rate at the fluid source origin followed by stabilization towards a constant rate. This might be the case of, for instance, breaking an initially sealed and highly-pressurized reservoir. Of course, after a certain period of injection at approximately constant rate, the rate of injection has to decrease until the pressure at the initially highly-pressurized fluid source equilibrates the fluid pressure of the surroundings; we neglect the rupture growth during that stage and also after injection ceases. Our results for sustained injection at constant volume rate can be thus discussed in the context of natural fluid-driven seismicity in seismic swarms episodes and aftershock sequences.
Indeed, the previous discussion on injection-induced aseismic slip in fractured rock masses is easily extendable to fault zones. Notably, seismicity is expected to be now constrained into a relatively welldefined fault plane of thickness that equals the size of the damage zone. Observed seismicity would be the result of instabilities that occur either in the main fault plane or in the fracture network of the damage zone. Similarly to the case of a fractured rock mass, the triggering of seismicity by aseismic-slip stress transfer would depend on how critically stressed the main fault plane or the damage zone fracture network is. The square root time migration of seismicity might be insufficient to discriminate whether aseismic slip or elevated fluid pressure is the direct triggering mechanism. Likewise, estimates of fault hydraulic diffusivity α based on seismicity patterns might rather represent the quantity αλ 2 .
Other phenomena that are potentially related to the diffusion of fluid pressure are silent earthquakes or slow slip events. Their occurrence is spatially well-correlated to predominantly frictionally stable regions in subduction zones with highly-pressurized fluids that operate at almost lithostatic pressures (Kato et al. 2010, Frank et al. 2015. Metamorphic dehydration reactions are likely to occur in these regions and have been invoked as a possible driving mechanism for slow slip events (Kato et al. 2010). To distinguish between possibly different physical mechanisms underlying slow and regular earthquakes, seismologists often search for scaling relations between the scalar moment M 0 and the duration of the event D (Gomberg et al. 2016). Assuming that slow slip events are driven by a sustained fluid source characterized by a constant injection rate over a fault patch that is small in comparison to the rupture area, and the rupture is allowed to grow unboundly, we obtain, from Eq. (38), that M 0 ∝ D 3/2 . This scaling relation differs from other relations proposed for a diffusion-controlled slow slip event model that follows M 0 ∝ D (Ide et al. 2007).

Permeability variations
Our model assumes that fluid flow occurs within a frictional interface characterized by a constant hydraulic transmissivity. However, permeability changes due to variations of the effective normal stress or, equivalently, the normal interfacial deformation/closure, are well-documented in the fracture/joint rock mechanics literature (e.g., Bandis et al. 1983) and fault mechanics literature as well (e.g., Rice 1992). In addition, fracture/fault dilatant-behavior can also induce significant permeability variations (e.g., Ciardo & Lecampion 2019). The effect of such hydro-mechanical couplings on the propagation of aseismic slip remains to be investigated in 3-D and requires the solution of the fully-coupled hydro-mechanical problem as solved, for instance, by Ciardo & Lecampion (2019) in 2-D.

Summary and concluding remarks
We have studied the quasi-static propagation of aseismic fault slip driven by fluid pressure diffusion under two different injection scenarios, namely, at constant volume rate and at constant pressure. Our model considers a frictional shear crack that grows in mixed mode (II+III) on a planar fault interface that separates two identical half-spaces of a three-dimensional, isotropic, homogeneous, linear elastic and impermeable solid. The fault interface is characterized by: a shear strength that is equal to the product of a constant friction coefficient and the local effective normal stress, a uniform stress state before injection, and a uniform and constant hydraulic transmissivity. The problem admits analytical treatments for circular ruptures which occur in the limit of a Poisson's ratio ν = 0, and it is solved numerically for the more general case in which the frictionally-constrained crack shape is to be determined as part of the solution (ν = 0).
For injection at constant volume rate from a point source, the fault rupture is self-similar. For the limiting case of a circular crack (ν = 0), the rupture radius evolves simply as R(t) = λL(t), where L(t) = √ 4αt is the nominal position of the fluid pressure front and λ is an amplification factor which is similar to the one presented by Bhattacharya & Viesca (2019) in their 2-D model. We derived an analytical solution for λ as a function of a unique dimensionless parameter T . The fault stress parameter T varies between 0 and +∞ and contains the information related to the pre-injection fault stress state and the parameters of the injection protocol. Whenever λ > 1, the rupture front outpaces the fluid pressure front. As in previous studies (Garagash & Germanovich 2012, Bhattacharya & Viesca 2019, two end-member cases have been identified, namely, critically stressed faults (T → 0) that largely outpace the fluid pressure front (λ 1), and marginally pressurized faults (T → +∞) that significantly lags the fluid pressure front (λ 1). Simple closed-form asymptotic expressions have been derived for λ and also for the axisymmetric slip distribution δ(r, t), for the two end-member cases. Other results include the rupture speed that decays with the square root of time as V r (t) = λ √ α/ √ t, and the accumulated fault slip at the injection point which is δ(r = 0, t) ≈ 3.5 (f ∆p * /µ) L(t) for critically stressed faults, and δ(r = 0, t) = (8/π) (f ∆p * /µ) R(t) for marginally pressurized faults, where f is the friction coefficient, ∆p * is the characteristic overpressure of the injection, and µ is the shear modulus.
For the more general case in which the Poisson's ratio is different than zero, we solved the problem of determining the equilibrium shape of the frictional shear crack over the entire parametric space. We found that the crack shape is quasi-elliptical and the aspect ratio is upper and lower bounded by 1/(1−ν) and (3 − ν)/(3 − 2ν). The two bounds are associated with the limiting cases of critically stressed faults and marginally pressurized faults, respectively. There is thus a strong dependence of the aspect ratio not only on the Poisson's ratio but also on the initial stress state and the driving force itself. Moreover, we found that the rupture area is Poisson's ratio-independent and grows simply as A r (t) = 4παλ 2 t. If λ > 1, the rupture area is greater than the diffusively pressurized fault patch. Interestingly, λ is the same amplification factor that for the circular rupture case, meaning that knowing the solution of the circular shear crack is sufficient to determine the area of any other resulting crack shape for any value of the Poisson's ratio and the same value of T . In addition, simple closed-form asymptotic expressions are provided for the semi-major a(t) and semi-minor b(t) axes of the quasi-elliptical crack that fully define the rupture front for the corresponding end-member cases.
For injection at constant pressure from a finite source of radius r w , the fault rupture is not self-similar. The rupture radius grows at large times as R(t) = λ(t) √ c 1 αt, where c 1 ≈ 2.25, √ c 1 αt is the nominal position of the fluid pressure front and λ(t) is an amplification factor known as function of dimensionless time αt/r 2 w and T . The fault stress parameter T varies in this case between 0 and 1. λ(t) decreases monotonically with time and the rupture radius expands as R(t) ∝ (αt) (1−T )/2 . For critically stressed faults (T → 0 =⇒ (1 − T ) /2 → 1/2), the rupture evolves almost diffusively, whereas for marginally pressurized faults (T → 1 =⇒ (1 − T ) /2 → 0), the rupture grows extremely slow. Generally speaking, the rupture front propagates slower than the diffusive fluid pressure front. Yet in some cases the rupture front outpaces the fluid pressure front, the latter will eventually catch the former if injection is sustained for a sufficient time.
Among the two injection scenarios considered, injection at constant volume rate is the one with broader implications. This is due to injection protocols in the geo-energy industry normally consist of a series of injections at constant volume rate, whereas naturally injected fluids into the Earth's crust are likely represented by the same kind of source. Since aseismic ruptures expand diffusively (proportional to the square root of time) for that type of injection, irrespective of the pre-injection stress state and the parameters of the injection, current interpretations of fluid-driven seismicity might need to be revisited.
Indeed, it is commonly assumed that seismicity clouds are driven by the direct effect of fluid pressure increase whenever seismic events are observed to spread away from the injection zone with square root time dependence (Shapiro et al. 1997, Bosl & Nur 2002, Parotidis et al. 2005, Chen et al. 2012, Hainzl et al. 2016, Ross et al. 2017, 2020. Our results challenge that interpretation and open the possibility that those episodes might be controlled by the stress transfer of a propagating aseismic rupture instead (Bhattacharya & Viesca 2019). This would be notably the case of critically stressed fractures/faults in which the rupture front is predicted to be systematically ahead of the fluid pressure front (λ 1). Furthermore, current estimates of reservoir and fault zone hydraulic diffusivity α based on seismicity patterns (e.g., Shapiro et al. 1997, Ross et al. 2017) might be rather related to the quantity αλ 2 , with λ being a representative amplification factor of the fractured rock mass or fault zone.
Another important finding is related to the scalar moment M 0 due to purely aseismic (quasi-static) motion. We found that it scales to the injected volume of fluid V as M 0 ∝ V 3/2 . Interestingly, this relation is the same as the one found by Galis et al. (2017) for self-arrested injection-induced seismic (dynamic) ruptures. Similarly, the seismic moment scales to the duration of fluid-driven aseismic ruptures D, approximately as M 0 ∝ D 3/2 . This scaling relation might be useful to interpret seismological observations of slow slip events driven by fluid sources such as dehydration reactions or others related to fault valving behavior, and it differs from other diffusion-controlled slow slip event models that have been proposed to follow M 0 ∝ D (Ide et al. 2007).
We expect our analytical and numerical results to provide a conceptual and quantitative framework to undertand various applied problems in geomechanics and geophysics associated with aseismic fracture/fault slip induced by fluid motion. Moreover, our analytical results provide a simple mean for verifying and benchmarking 3-D numerical solvers as performed here.
CRediT authorship contribution statement

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Funding
The results were obtained within the EMOD project (Engineering model for hydraulic stimulation). The EMOD project benefits from a grant (research contract SI/502081-01) and an exploration subsidy (contract number MF-021-GEO-ERK) of the Swiss federal office of energy for the EGS geothermal project in Haute-Sorne, canton of Jura, which is gratefully acknowledged. A. Sáez was partially funded by the Federal Commission for Scholarships for Foreign Students via the Swiss Government Excellence Scholarship. P. Bhattacharya was supported by a start-up research grant from the National Institute of Science Education and Research, Bhubaneswar. R. C. Viesca was supported by NSF award 1653382.

A The consistent tangent operator
The Newton-Raphson iterations of the backward Euler time integration scheme require the computation of the corresponding Jacobian matrix J ∈ R N ×N . By differentiating the residual form of Eq. (9), the Jacobian matrix is simply given by J = E H +C T O , where C T O = −∂∆t /∂∆d is the so-called consistent tangent operator of elastoplasticity. In order to derive an analytical expression for C T O , we consider the consistency condition ∆γḞ = 0, which states that when plastic flow occurs (i.e., ∆γ > 0), the stress state t has to remain on the yield function and thusḞ = 0. This additional equation can be read in incremental form as
Note that the consistent tangent operator C T O is a block diagonal matrix and is composed by blocks that we denote as C m T O ∈ R 3×3 . There is one C m T O matrix for every m-th collocation point. By combining Eqs. (10) to (13) plus the previous consistency condition, we derive the following expression for C m where θ = arctan (t 2 /t 1 ) with t 1 and t 2 the two local components of the shear traction vector at the m-th collocation point (the superscript m is omitted), and k 1 , k 2 and k 3 are the corresponding entries of the diagonal elastic stiffness matrix D.
Note that C m T O is a null matrix if ∆γ m = 0 (i.e., if the collocation point state is elastic or, in other words, no slip has occured).

B Propagation condition for a constant-friction circular shear crack under axisymmetric shear load
The stress intensity factors of a circular crack of radius R under an arbitrary shear traction vector of components σ xz and σ yz (with regard to the reference frame showed in Fig. 1) applied anti-symmetrically on the crack surfaces read as (Fabrikant 1989, Lai et al. 2002) where φ is the polar angular coordinate, such that tan (φ) = x/y.
Consider a shear load of axisymmetric magnitude ∆τ (r) along the x direction, such that σ xz (r, θ) = ∆τ (r) and σ yz (r, θ) = 0. Evaluating the integral of the right-hand side of the previous equation with regard to θ, we obtain: Let us consider the energy release rate of a pure shear crack in 3-D elasticity, G = K 2 II /E + K 2 III /2µ (Lawn 1993), where E is the plane strain modulus. Using the previous equations for the stress intensity factors in the limiting case of a material with Poisson's ratio ν = 0 (E = E = 2µ), we obtain the following expression for the energy release rate, The fracture energy G c of a constant-friction shear crack is zero, such that Grifith energy balance G = G c reduces simply to R 0 ∆τ (r) √ R 2 − r 2 rdr = 0, that is the expression used in the main text as the condition for crack propagation.

C Asymptotics of fault slip for circular ruptures driven by injection at constant volume rate
The quasi-static elastic equilibrium that relates the fault slip distribution δ to the shear stress drop ∆τ within an axisymmetric circular shear crack (ν = 0) of radius R(t) is (Salamon & Dundurs 1971, 1977 ∆τ (r, t) = µ 2π where K and E are the complete elliptic integrals of the first and second kind, respectively, and k(x) = 2 √ x/ (1 + x).
The inverse relation of the previous equation is given by Sneddon (1951) wherer = r/R(t) is the self-similar radial coordinate.
Considering that the shear stress drop ∆τ (r, t) = τ 0 −f σ o − ∆p * E 1 r 2 /4αt (Eq. (19)), the fault stress parameter T = (f σ o − τ o ) /f ∆p * (Eq. (17)), and the amplification factor is defined as λ = R(t)/L(t) with L(t) = √ 4αt, we can recast the above equation in dimensionless form, whereδ (r; T ) is the normalized self-similar slip distribution that is unique for a given value of T . We recall that the amplification factor λ is known by Eq. (21) as a function of T as well.
The latter is indeed valid for r L(t) only. It corresponds to the "outer" solution of the critically stressed limit in which the reduction of frictional strength due to the fluid pressure perturbation can be approximated as a point force. Eqs. (53) and (54) can be equivalently derived by using the asymptotic expressions of the fluid pressure perturbation in the limiting cases (see details on those approximations in section 3.2). In addition, the condition for having no singularity at the crack tip, ∂δ/∂r = 0 at r = R, is equivalent to the crack propagation condition, Eq. (18), and will lead to same expressions that relates λ and T , Eqs. (22) and (23), in the corresponding end-member cases.
Eqs. (53) and (54) are plotted in Fig. 14 together with the slip profiles obtained from the numerical simulations for values of T that are representative of the limiting cases. We use logarithmic scale in the critically stressed limit in order to facilitate the comparison. Note that the marginally pressurized limit is reached for values of T 4 (Fig 14b). On the other hand, in the critically stressed limit (Fig 14a), the "outer" solution breaks for small r (it diverges at r = 0 indeed); an "inner" solution should be derived in order to properly approximate the slip distribution in the domain in which r ∼ O (L(t)).

D Analytical solution in 2-D for a frictional shear crack driven by injection at constant volume rate
Consider the same problem formulated in section 2.1 but in 2-D elasticity. Fluid is injected via a point source at x = 0 into a planar 1-D frictional interface located along the x-axis. Injection is sustained for t > 0 at a constant volume rate q w [L/T ] per unit fault thickness and unit out-of-the-plane length. The quasi-static crack propagation condition for a 1-D straight constant-friction (zero-toughness) shear crack (either mode II or III) of half-crack length (t) reduces to (Barenblatt 1962) (t) where ∆τ (x, t) is the shear stress drop given by ∆τ (x, t) = τ 0 − f (σ o − ∆p(t)Π(ξ)) , with ξ = x/ √ 4αt and ∆p(t) = 2q w η √ πk √ αt, Π(ξ) = exp −ξ 2 − √ π |ξ| Erfc (|ξ|) .
As defined in the main text, η is the fluid dynamic viscosity, k is the fault intrinsic permeability, and α is the fault hydraulic diffusivity.
Let us define the nominal position of the fluid pressure front d (t) = √ 4αt, a time-dependent amplification factor λ in the form λ(t) = (t)/ d (t), and the following characteristic length

The crack propagation condition can be then rewritten in dimensionless form
that represents a unique relation between the amplification factor λ and the ratio between the position of the fluid pressure front d and the characteristic length x c .
The left-hand side of the previous integral can be evaluated analytically to obtain the following implicit equation for λ as a function of d /x c , exp −λ 2 /2 1 + λ 2 I 0 λ 2 /2 + λ 2 I 1 λ 2 /2 − 2λ/ √ π = x c / d , where I 0 and I 1 are the Bessel functions of the first kind of zero and first order, respectively.
The previous equation is the one reproduced in the main text and it is plotted in Fig. 13. Asymptotic expansions of the left-hand side of this equation for small and large λ lead to the following closed-form asymptotic solutions for the short-run-out rupture (λ 1) and long-run-out rupture (λ 1) regimes: , for λ 1 , that are also plotted in Fig. 13.
From the asymptotic for the short-run-out rupture regime, we note that solutions are defined only for d /x c ≥ 1. The limit of d /x c = 1 represents the instant in which activation of slip (or crack nucleation) happens and implies that d = x c . From the latter, it becomes clear that the characteristic length x c is the size of the pressurized patch necessary to the crack starts growing.
On the other hand, from the asymptotic for the long-run-out rupture regime and the definition of λ, we can write Hence, the ultimate behavior of the crack is a steady propagation regime at constant rupture speed V R equal to which is the other equation used in the main text.