Hostname: page-component-7b7567c8c-7pwfn Total loading time: 0 Render date: 2026-11-12T18:24:17.337Z Has data issue: false hasContentIssue false
  • English
  • Français

Asymptotic solutions for self-similarly expanding fault slip induced by fluid injection at constant rate

Published online by Cambridge University Press:  17 September 2025

Robert C. ViescaOpen the ORCID record for Robert C. Viesca[Opens in a new window]
Robert C. Viesca*
Affiliation:
Department of Civil and Environmental Engineering, Tufts University, Medford, MA 02155, USA
*
Corresponding author: Robert C. Viesca, robert.viesca@tufts.edu
Rights & Permissions [Opens in a new window]

Abstract

We examine the circular, self-similar expansion of frictional rupture due to fluid injected at a constant rate. Fluid migrates within a thin permeable layer parallel to and containing the fault plane. When the Poisson ratio $\nu =0$, self-similarity of the fluid pressure implies fault slip also evolves in an axisymmetric, self-similar manner, reducing the three-dimensional problem for the evolution of fault slip to a single self-similar dimension. The rupture radius grows as $\lambda \sqrt {4\alpha _{hy} t}$, where $t$ is time since the start of injection and $\alpha _{hy}$ is the hydraulic diffusivity of the pore fluid pressure. The prefactor $\lambda$ is determined by a single parameter, $T$, which depends on the pre-injection stress state and injection conditions. The prefactor has the range $0\lt \lambda \lt \infty$, the lower and upper limits of which correspond to marginal pressurisation of the fault and critically stressed conditions, in which the fault-resolved shear stress is close to the pre-injection fault strength. In both limits, we derive solutions for slip by perturbation expansion, to arbitrary order. In the marginally pressurised limit ($\lambda \rightarrow 0$), the perturbation is regular and the series expansion is convergent. For the critically stressed limit ($\lambda \rightarrow \infty$), the perturbation is singular, contains a boundary layer and an outer solution, and the series is divergent. In this case, we provide a composite solution with uniform convergence over the entire rupture using a matched asymptotic expansion. We provide error estimates of the asymptotic expansions in both limits and demonstrate optimal truncation of the singular perturbation in the critically stressed limit.

JFM classification

Boundary Layers: Boundary LayersGeophysical and Geological Flows: Geophysical and Geological Flows

Information

Type
JFM Papers
Information
Journal of Fluid Mechanics , Volume 1019 , 25 September 2025, A27
Check for updates
Copyright
© The Author(s), 2025. Published by Cambridge University Press

1. Introduction

Boundary layers and matched asymptotic perturbations are implicit in many problems in fracture mechanics, but details are often subsumed within the so-called stress intensity factor. This factor provides the leading-order matching condition between an outer solution on the scale of the entire fracture and an inner solution representing a physical – often weakening – process in the vicinity of the crack tip. When the small-scale-yielding condition for the separation of scales between the entire fracture and this near-tip cohesive zone is met, the details of the breakdown process are often neglected, yet assumed to be universal for a given material such that a critical value of the stress intensity factor for a propagating crack is assumed to be a material constant. As a consequence, a cohesive-zone boundary layer solution near the tip is seldom explicitly sought and even more rarely matched with an outer solution to provide a composite solution of crack-face displacement along the entire crack. In practice, this step is often unnecessary, but may be a missed opportunity to demonstrate a more fundamental understanding of models involving fracture and to exhibit perturbation methods as applied to fracture.

There nonetheless has been a substantial series of works that have pursued singular perturbations involving boundary layers in classical fracture problems and in those coupling fracture with other physical processes. In a series of works, Edmunds & Willis (Reference Edmunds and Willis1976a , Reference Edmunds and Willisb , Reference Edmunds and Willis1977) examined higher-order matching of the far-field with the near-field of a crack-tip cohesive zone. Among other results, the authors demonstrated a correction to the crack length appearing in standard expressions for the stress-intensity factor, the corrections being equivalent to those of second-order to the energy release rate. For the specific geometry of a long, thin crack running parallel to a free surface, as may occur in spalling, several workers examined an inner solution corresponding to a semi-infinite crack as well as its matching with an outer solution corresponding to the thin ligament deforming as a simple Euler–Bernoulli beam (Zlatin & Krapkov Reference Zlatin and Khrapkov1986; Thouless et al. Reference Thouless, Evans, Ashby and Hutchinson1987; Dyskin, Germanovich & Ustinov Reference Dyskin, Germanovich and Ustinov2000). Lister, Skinner & Large (Reference Lister, Skinner and Large2019) performed similar analysis considering the more elaborate case where the crack is being driven by the migration of a viscous fluid. Such hydraulic fracture in its own right has been a source of novel problems in the asymptotic, boundary-layer analysis of fractures. Early work included self-similar solutions and asymptotic examination of near-tip viscous flow (Spence & Sharp Reference Spence and Sharp1985; Lister Reference Lister1990; Desroches et al. Reference Desroches, Detournay, Lenoach, Papanastasiou, Pearson, Thiercelin and Cheng1994; Garagash & Detournay Reference Garagash and Detournay2000). This was followed by matched-asymptotic expansions of self-similarly expanding hydraulic fractures in the so-called viscosity-dominated regime, in which, to leading order, effects of fracture toughness are confined to a near-tip boundary layer (Garagash & Detournay Reference Garagash and Detournay2005), as well as regular asymptotic expansions in the so-called large-toughness regime (Savitski & Detournay Reference Savitski and Detournay2002; Garagash Reference Garagash2006). Boundary-layer solutions also appear in problems of the fracture of poroelastic media. Atkinson & Craster (Reference Atkinson and Craster1994) examined the sudden appearance of a finite, stress-free crack in a porous medium. This problem has an early-time solution that can be considered to be the composition of an ‘undrained’ outer-solution over most of the fracture, in which fluid flow is neglected, as well as an inner solution valid near the crack tips, the surroundings of which are in an effectively ‘drained’ state and for which fluid flow is sufficiently quick as to instantaneously establish a constant, uniform fluid pressure. An earlier series of works also considered the steady propagation of a semi-infinite fracture in a poroelastic medium for a variety of permeability conditions along the fault, and in which boundary layers appear in the vicinity of the crack tip (Simons & Rice Reference Simons and Rice1976; Simons Reference Simons1977; Rudnicki Reference Rudnicki1990; Rudnicki & Koutsibelas Reference Rudnicki and Koutsibelas1991; Craster & Atkinson Reference Craster and Atkinson1992,Reference Craster and Atkinson1994). In all of the above-mentioned works, perturbations are carried out to within two terms, with the exception of the work of Garagash (Reference Garagash2006) who provides results of numerical solutions for finitely higher-order expansions.

Here, we look to a model problem that provides the opportunity to examine both regular and singular perturbation methods to arbitrary order and to demonstrate optimal truncation of the latter. The model problem consists of a singular integral equation for a circular crack in a three-dimensional medium, the growth of which is driven by fluid injection at constant volume rate into a narrow permeable region surrounding a planar geologic fault, such that fluids permeate away from the injection site and remain within the fault plane. The injection elevates pore fluid pressure along the fault and reduces the effective normal stress carried by solid contacts, and hence the frictional strength. We will investigate an elementary scenario in which the initial stress state is uniform and the fault’s frictional strength follows a Coulomb friction law. Injection of fluid into the fault plane – and its subsequent migration along that plane – elevates fluid pressure there, which drives the loss of strength and permits the nucleation and growth of a crack. Because the functional form of the forcing is known a priori, we can directly express a self-similar solution for slip and crack radius in integral form using a solution of Sneddon (Reference Sneddon1951). This solution cannot be calculated directly and must be evaluated numerically. However, asymptotic solutions can be found using perturbation techniques. The problem has one parameter with two end-member limits: in one limit, we derive a regular perturbation expansion, while in the other, we derive a singular perturbation expansion accounting for the existence of a boundary layer. The boundary layer results from a separation of the length scale characteristic of the forcing driving the crack from the size of the crack itself. In the limit in which the boundary layer exists, we will not need to explicitly match the asymptotic expansions to connect the so-called inner solution of the boundary layer with the outer solution away from the layer. Rather, the inner and outer solutions we calculate will inherently satisfy matching conditions, owing to the use of the inversion formula. Previously, we performed a comparable analysis for the self-similar growth of a frictional mode-II or -III fracture with a constant friction coefficient under fluid injection at constant pressure, providing the first two terms in the perturbation series solution in each end-member limit (Viesca Reference Viesca2021); leading-order asymptotics in these limits were also sought and solved for by Garagash & Germanovich (Reference Garagash and Germanovich2012) via a mix of closed-form and numerical techniques for the case of a mode-II or -III fracture under the same injection conditions, except with slip-weakening friction, which precluded self-similarity but permitted dynamic rupture nucleation. Here, we address the mixed-mode problem of a circular shear crack.

This problem was previously considered numerically by Sáez et al. (Reference Sáez, Lecampion, Bhattacharya and Viesca2022) who considered the more general free-boundary problem in which the crack boundary is not known a priori and part of the solution. When the Poisson ratio $\nu =0$ , the crack boundary is expected to be circular and grow in a self-similar manner. In this axisymmetric case, numerical solutions were provided spanning the range of the problem parameter. In addition, leading-order asymptotic solutions in this circular-crack case were provided in the limit of two end-member values of the problem parameter: a complete leading-order solution was provided in one limit, and an outer solution provided in the other. A boundary element method was employed to solve the free-boundary problem when $\nu \gt 0$ and the numerical axisymmetric solutions, in addition to the asymptotic behaviour, were used to verify boundary-element solutions in the $\nu =0$ limit. The solutions of Sáez et al. (Reference Sáez, Lecampion, Bhattacharya and Viesca2022) demonstrated that, for cases in which $\nu \neq 0$ , crack growth also proceeds self-similarly and its boundary is nearly elliptical. These solutions also demonstrated that the growth of the elliptical rupture area ( $\nu \neq 0$ ) was, to within the accuracy of the solutions, equal to that of the circular crack ( $\nu =0$ ) for a given value of the problem parameter. This has led to a renewed interest in using a circular-crack approximation as a reduced-order model for related problems (Sáez & Lecampion Reference Sáez and Lecampion2023, Reference Sáez and Lecampion2024; Sáez, Passeègue, & Lecampion, Reference Sáez, Passelègue and Lecampion2025), as well as a benchmark. Here, we return to this problem and provide the complete asymptotic expansions in both limits of the problem parameter, including the full inner and outer solutions in the limit for which singular perturbations occur and a boundary layer exists.

Independent of its perturbation aspects, this problem is of interest to understand the activation of geologic faults in the subsurface in response to transiently elevated pore fluid pressure. The fundamental and industrial interest for such problems range from understanding the interplay between natural fluid migration and aseismic or seismic fault activity to monitoring of the artificial injection of fluids, such as to avoid the activation of faults to preserve the integrity of an imperable caprock formation, or to intentionally enhance the permeability of naturally occurring faults for the extraction of geothermal energy by inducing their slip. In any context, non-numerical solutions for models of such problems are rare. Such solutions provide fundamental scaling relations for the propagation of the rupture front and the accumulation of fault slip, in addition to its distribution, as they may vary with the strength of the injection and the pre-existing state of stress on the fault. Comprehensive analytical solutions are typically only available when model complexity is low, as is the case for the elementary model examined here. Nonetheless, these solutions provide a basis to understand how increasing model complexity alters model behaviour. Areas of potential complexity in these models may include representations of the fault hydrology, friction or geometry. For example, several prior works have examined the evolution of fault slip following fluid injection assuming friction laws alternative to the Coulomb law assumed here, such as slip-weakening or rate- and state-dependent friction (e.g. Garagash & Germanovich Reference Garagash and Germanovich2012; Viesca & Rice Reference Viesca and Rice2012; Dublanchet Reference Dublanchet2019; Im & Avouac Reference Im and Avouac2024). These have predominately examined ruptures within a two-dimensional medium, such as for in-plane or anti-plane deformation.

The outline of this work is as follows. In § 2, we formulate the problem determining the fault slip distribution and rupture front location. We provide an integral representation of the solution for these quantities. We provide asymptotic expressions for the solution of the rupture front location as it depends on the pre-injection state of fault traction and strength and a diffusive length scale $\sqrt {\alpha t}$ . In § 3, we provide singular perturbation expansions of this solution in the so-called critically stressed limit of the problem parameter. In this limit, the solution for slip contains a boundary layer at the centre of the crack, about the injection point, over distances comparable to the diffusive length scale. We provide the expansion for this inner solution to arbitrary order in the perturbation parameter. Over distances comparable to the crack radius but away from the diffusive boundary layer, the details of fluid diffusion are negligible, and we likewise provide the expansion for the so-called outer behaviour of the solution to arbitrary order. This section demonstrates the asymptotic matching of the inner and outer solutions, calculates their overlap, and constructions a composite solution that uniformly approximates the solution over the entire rupture. In § 4, we provide regular perturbation expansion of the solution for slip in the so-called marginally pressurised limit of the problem parameter. In § 5, we examine the accumulation of the peak slip and aseismic moment across all values of the problem parameter and provide asymptotic expansions for both quantities in the critically stressed and marginally pressurised limits. In § 6, we provide error estimates of the singular and regular perturbation series in the two limits and demonstrate optimal truncation of the singular perturbation series. Finally, in § 7, we discuss several issues concerning the presented solutions, including the choice of perturbation parameter, the importance of the solution to a problem for which $\nu =0$ and relevant applications.

2. Problem formulation

For a circular shear fracture of radius $R(t)$ undergoing relative slip $\delta (r,t)$ in a single direction, the shear traction on the crack faces $\tau (r,t)$ , opposing the direction of slip, is given by (Salamon & Dundurs Reference Salamon and Dundurs1971,Reference Salamon and Dundurs1977; Bhattacharya & Viesca Reference Bhattacharya and Viesca2019)

(2.1) \begin{equation} \tau (r,t)=\tau _b+\frac {\mu }{2\pi }\int _0^{R(t)} \left (\frac {E[k(r/s)]}{s-r}+\frac {K[k(r/s)]}{r+s}\right )\frac {\partial \delta (s,t)}{\partial s}\,{\rm d}s, \end{equation}

where $k(u)=2\sqrt {u}/(1+u)$ , and $E$ and $K$ are the complete elliptical integrals of the second and first kind. The background shear stress resolved on the plane of the fracture in the absence of slip, $\tau _b$ , is here taken to be uniform. The relation (2.1) presumes that the Poisson ratio $\nu =0$ . This expression is equivalent to that relating relative crack opening and crack-face tractions for a circular mode-I crack. In that case, this expression holds for any value of the Poisson ratio, provided that $\mu$ here is replaced by $\mu /(1-\nu )$ . The inversion of this equation for the relative displacement $\delta$ was done in the context of mode-I equation (Sneddon Reference Sneddon1951), which we adapt here for our particular mixed-mode II/III case with $\nu =0$ ,

(2.2) \begin{equation} \delta (\rho ,t)=\frac {4}{\pi }\frac {R(t)}{\mu }\int _\rho ^1\frac {\eta }{\sqrt {\eta ^2-\rho ^2}}\int _0^1\frac {x\,\Delta \tau (x\eta R,t)}{\sqrt {1-x^2}}\,{\rm d}x\,{\rm d}\eta, \end{equation}

where

(2.3) \begin{equation} \Delta \tau (r,t)=\tau _b-\tau (r,t) \end{equation}

and $\rho =r/R(t)$ .

We may rearrange this inversion to isolate the contribution of a crack-tip singularity to the slip distribution. To do so, we first denote the innermost integral

(2.4) \begin{equation} g(\eta )=\int _0^1\frac {x\,\Delta \tau (x\eta R,t)}{\sqrt {1-x^2}}\,{\rm d}x, \end{equation}

where for brevity and without loss of generality, we have suppressed explicit mention of the time dependence of $g$ . The inversion (2.2) reduces to

(2.5) \begin{equation} \delta (\rho ,t)=\frac {4}{\pi }\frac {R(t)}{\mu }\int _\rho ^1\frac {\eta }{\sqrt {\eta ^2-\rho ^2}}g(\eta )\,{\rm d}\eta .\end{equation}

Integration of (2.5) by parts yields

(2.6) \begin{equation} \frac {\delta (\rho ,t)}{4 R(t)/(\pi \mu )}=g(1)\sqrt {1-\rho ^2}-\int _\rho ^1\sqrt {\eta ^2-\rho ^2}\,\frac {{\rm d}g}{{\rm d}\eta }\,{\rm d}\eta, \end{equation}

where we identify that

(2.7) \begin{equation} g(1)=\int _0^1 \frac {x \Delta \tau (x R,t)}{\sqrt {1-x^2}}\,{\rm d}x \end{equation}

is a stress intensity factor as it is a pre-factor for the relative-displacement field behind the crack tip, $\delta \sim g(1)\sqrt {1-r}$ , corresponding to a singular stress field ahead of the crack tip, $\tau \sim g(1)/\sqrt {r-1}$ . This shear stress intensity factor is identically zero for problems in which the shear stress must remain bounded, without singularity, owing to the finite strength of the fault, the source for which is Coulomb friction.

We now look to consider a particular form of $\Delta \tau$ corresponding to the injection of fluid at a constant volumetric rate $Q$ into a thin fault zone whose Coulomb friction coefficient remains at a constant value, $f$ . The strength of the fault is given by the product of $f$ and the effective normal stress $\sigma -p$ ,

(2.8) \begin{equation} \tau _s=f(\sigma -p(r,t)), \end{equation}

where $\sigma$ is the total normal stress, which is constant and uniform along the fault, $p$ is the local fluid pressure at a point on the fault, which is presumed to have an initial, uniform value $p_o$ plus a perturbation due to injection, given by the solution to the axisymmetric diffusion equation governing fluid pressure in the plane of the thin fault zone

(2.9) \begin{equation} \frac {\partial p}{\partial t}=\alpha _{hy}\frac {1}{r}\frac {\partial }{\partial r}\left (r\frac {\partial p}{\partial r} \right ), \end{equation}

where $\alpha _{hy}=k/(\eta _f \beta )$ is the hydraulic diffusivity, $k$ is the fault-zone permeability, $\eta _f$ is the fluid viscosity, and $\beta$ is the composite compressibility of the porous matrix and fluid. The corresponding solution of interest is that subject to the injection condition

(2.10) \begin{equation} -2\pi r k\frac {\partial p}{\partial r}=\frac {Q}{h} \end{equation}

and is given by

(2.11) \begin{equation} p(r,t)=p_o+\Delta p E_1 \big[r^2/(\alpha t) \big], \end{equation}

where $E_1$ is the exponential integral, $\alpha = 4\alpha _{hy}$ and $\Delta p = Q/(4\pi k h)$ . The thinness of the fault zone is relevant when considering that the fault zone is a permeable, poroelastic conduit confined within an elastic medium, which may have a stiffness dissimilar to the fault zone. Here, we consider the case where the fault-zone thickness $h$ is much smaller than the diffusive length scale $\sqrt {\alpha t}$ . Under this condition, there is no change in the fault-normal stress over distances comparable to $\sqrt {\alpha t}$ due to localised, constrained swelling of the fault zone around the injection source. In this case, fluid pressure is decoupled from volumetric straining of the fault zone and follows the autonomous diffusion equation (Marck, Savitski & Detournay (Reference Marck, Savitski and Detournay2015) consider the injection scenario here; see also Jacquey & Viesca (Reference Jacquey and Viesca2023), for similar results under plane-strain conditions; earlier works, such as Garagash & Germanovich (Reference Garagash and Germanovich2012), have implicitly operated in this limit).

Self-similarity of the pore-fluid pressure suggests self-similarity of the rupture front and slip profile. In this case, the rupture radius may be written as

(2.12) \begin{equation} R(t)=\lambda \sqrt {\alpha t}, \end{equation}

where $\lambda$ is a pre-factor to be determined. The non-singular stress condition, $f(1)=0$ , provides a condition to solve for the pre-factor $\lambda$ in terms of other problem parameters. This condition is imposed because the finite-frictional strength of the fault precludes a stress singularity. To impose this non-singular condition, we note that the frictional strength condition requires fault shear stress to equal fault shear strength, $\tau (r,t)=\tau _s(r,t)$ , when and where sliding occurs. Combining this equation with (2.3), and (2.7)–(2.12), and rearranging,

(2.13) \begin{equation} T=\int _0^1 \frac {x\, E_1\!\left [(\lambda x)^2\right ]}{\sqrt {1-x^2}}\,{\rm d}x, \end{equation}

where the dimensionless parameter

(2.14) \begin{equation} T=\frac {\sigma '}{\Delta p}\left (1-\frac {\tau _b}{\tau _p}\right ) \end{equation}

is a measure of distance of the fault stress state from critically stressed conditions, $\tau _b\rightarrow \tau _p$ , in which $T\rightarrow 0$ and the pre-injection strength of the fault is $\tau _p=f\sigma '$ , and also where $\sigma '=\sigma -p_o$ is the pre-injection effective fault normal stress. When the initial stress $\tau _b$ is a finite distance away from $\tau _p$ , but the injection-induced pressure coefficient $\Delta p$ is small compared with $\sigma '$ , the fault is marginally pressurised and $T\rightarrow \infty$ .

We may anticipate the behaviour of $\lambda$ in these two end-member regimes: i.e. we expect that the rupture will progress quickly relatively to pore-fluid pressure diffusion ( $\lambda \gg 1$ ) when the fault is critically stressed and relatively slowly ( $\lambda \ll 1$ ) when marginally pressurised. While a full solution to (2.13) was stated by Sáez et al. (Reference Sáez, Lecampion, Bhattacharya and Viesca2022) as an infinite (hypergeometric) series, we can derive concise expressions of the relation between $T$ and $\lambda$ in these limits by asymptotically expanding the integral (2.13).

For large values of $\lambda$ , we apply the change of variable $u=x\lambda$ ,

(2.15) \begin{equation} T=\frac {1}{\lambda ^2}\int _0^\lambda \frac {u\, E_1\!\left [u^2\right ]}{\sqrt {1-(u/\lambda )^2}}\,{\rm d}u ,\end{equation}

followed by the Taylor expansion

(2.16) \begin{equation} \frac {1}{\sqrt {1-(u/\lambda )^2}}=1+\frac {(u/\lambda )^2}{2}+\cdots +\frac {1}{\sqrt {\pi }}\frac {[(n-1)/2]!}{(n/2)!}(u/\lambda )^n \end{equation}

for even $n$ , and

(2.17) \begin{equation} T= \sum _{n=0,2, \ldots } ^\infty \frac {1}{\lambda ^{2+n}} \frac {1}{\sqrt {\pi }}\frac {[(n-1)/2]!}{(n/2)!} \int _0^\lambda u^{n+1}\, E_1 \big(u^2 \big){\rm d}u. \end{equation}

Passing the limit of the integrals in (2.17) to infinity without loss of generality of the power-law asymptotic expansion, owing to the exponentially fast decay of the integrand, we may first evaluate the integral by parts, second by subsequently using another change of variable $w=u^2$ and third by using the definition of the gamma function $\varGamma (z)= \int _0^\infty s^{z-1}e^{-s} \,{\rm d}s$ ,

(2.18) \begin{align} \begin{split} \int _0^\infty u^{n+1}\, E_1\!\left (u^2\right ){\rm d}u &= \frac {2}{n+2}\int _0^\infty u^{n+1} e^{-u^2} {\rm d}u \\[3pt] &=\frac {1}{n+2}\int _0^\infty w^{n/2} e^{-w} {\rm d}w=\frac {\varGamma (n/2+1)}{n+2}=\frac {(n/2)!}{n+2}. \end{split} \end{align}

Combining (2.17) with (2.18), and changing the index of summation $n=2m$ , we find the asymptotic expansion of $T$ for large $\lambda$ in the critically stressed limit,

(2.19) \begin{equation} T\sim \sum _{m=0}^\infty \frac {1}{\sqrt {\pi }}\frac {[(2m-1)/2]!}{2m+2} \frac {1}{\lambda ^{2+2m}}=\frac {1}{2\lambda ^2}+\frac {1}{8\lambda ^4}+\frac {1}{8\lambda ^6}+\cdots \end{equation}

For small values of $\lambda$ , we may evaluate (2.15) using the series expansion of the exponential integral for small values of its argument (e.g. Bender & Orszag Reference Bender and Orszag1999)

(2.20) \begin{equation} E_1 \big [(\lambda x)^2\big ]= -\gamma - 2 \ln (\lambda x)-\sum _{n=1}^\infty (-1)^n \frac {(\lambda x)^{2n}}{n \boldsymbol{\cdot }n!}, \end{equation}

where $\gamma =0.57721\ldots$ is the Euler–Maraschoni constant, such that

(2.21) \begin{align} T&= - \gamma \int _0^1\frac {x}{\sqrt {1-x^2}}\,{\rm d}x -2\int _0^1 \frac {\ln (\lambda x)}{\sqrt {1-x^2}}\,{\rm d}x-\sum _{n=1}^\infty (-1)^n \frac {\lambda ^{2n}}{n \boldsymbol{\cdot }n!}\int _0^1 \frac {x^{2n +1}}{\sqrt {1-x^2}}\,{\rm d}x \nonumber\\[3 pt] &= - \gamma + 2 -2 \ln (2\lambda ) - \sum _{n=1}^\infty (-1)^n \frac {\lambda ^{2n}}{n \boldsymbol{\cdot }n!}\left ( n! \frac {\sqrt {\pi }/2}{(n+1/2)!}\right ) \nonumber\\[3 pt] &= - \gamma + 2 -2 \ln (2\lambda ) - \sum _{n=1}^\infty (-1)^n \frac {\sqrt {\pi }/2}{n\boldsymbol{\cdot }(n+1/2)!} \lambda ^{2n} \nonumber\\[3 pt] &= - \gamma + 2 -2 \ln (2\lambda ) + \frac {2}{3}\lambda ^2-\frac {2}{15}\lambda ^4+\frac {8}{315}\lambda ^6 -\ldots \end{align}

in which the first integral is readily evaluated, the second likewise but by parts, and the evaluation of the last integral is to be found as an ancillary exercise in Appendix D, where we identify it as $I_{2n+1}$ . In figure 1, we compare the variation of $T$ with $\lambda$ , satisfying (2.15), against the marginally pressurised and critically stressed expansions.

We are interested in deriving series expansions for slip in the limits of large and small $\lambda$ . In the marginally pressurised limit, $\lambda \ll 1$ , the rupture radius lags far behind the fluid diffusive length scale $\sqrt {\alpha t}$ . In the critically stressed limit, where $\lambda \gg 1$ , the rupture radius is much larger than the diffusive length scale, $R(t)\gg \sqrt {\alpha t}$ , and with this scale separation, there will be an outer solution valid over distances ${\sim} R$ and an inner solution relevant over distances ${\sim} \sqrt {\alpha t}$ near $r=0$ . This inner solution is, in effect, an interior boundary layer. In what follows, we present in § 3 the asymptotic expansion of the inner and outer solutions for the critically stressed case, as well as their asymptotic matching and the construction of a composite solution; and we show in § 4 the asymptotic expansion of slip in the marginally pressurised limit. Details of the solution procedure will be relegated to Appendix A for the inner solution, Appendices B and C for the outer solution, and Appendix D for the marginally pressurised solution.

Figure 1.

Plot of the variation of the stress-injection parameter $T$ (solid black), which characterises the pre-injection state of fault stress and strength as well as the strength of the injection, with the prefactor $\lambda$ for the self-similar expansion of the crack radius, $R(t)=\lambda \sqrt {\alpha t}$ . The marginally pressurised (2.21, blue-dashed) and critically stressed (2.19, red-dashed) are the perturbation expansions in the small and large $\lambda$ limit, respectively. The marginally pressurised expansion is carried to $O(\lambda ^8)$ and the critically stressed expansion is carried to $O(\lambda ^{-4})$ , which represents an optimal truncation for $\lambda$ near 1 (discussed further in § 6.1).

3. Critically stressed faults: inner, outer and composite solutions

For small values of the stress parameter $T$ (large $\lambda$ ), the critically stressed fault is highly sensitive to the injection of fluids: rupture outpaces the migration of injected fluids. We perform a singular perturbation expansion of the self-similar slip profile in terms of the small parameter $\lambda ^{-1}$ . In the following, we provide the behaviour of slip on an interior boundary layer about the origin in § 3.1. This solution captures behaviour near the origin and will asymptotically decay away from it. Subsequently, in § 3.2, we provide an outer solution that captures behaviour at distances greater than the diffusive length away from the origin and up to the rupture front. The outer solution’s approximation breaks down at the origin, where that solution is singular. We will find that the manner of the asymptotic decay of the inner solution matches the singularity of the outer solution and that the two overlap over distances intermediate to the diffusive length scale and rupture front position. In § 3.3, we construct a composite solution as the sum of the inner and outer solutions, less the common overlap.

3.1. Inner solution

To examine the variation of slip over distances of the order of the diffusive length scale $\sqrt {\alpha t}$ away from the origin, we use the scaled radial position $s=r/\sqrt {\alpha t}$ , which has the domain $0\leqslant s \leqslant \lambda$ . While the details of the solution are deferred to Appendix A, the final result is that the asymptotic expansion for slip in the large $\lambda$ limit is given by $\delta \sim \delta _{in}$ , where

(3.1) \begin{align} \frac {\delta _{in}(s,t)}{\sqrt {\alpha t} f \Delta p/\mu }=\delta _0(s,t)+\frac {4}{\pi ^2}\sum _{j=0}^\infty \frac {1}{\lambda ^{2j+1}}\frac {1}{2j+1}\left [\sum _{k=0}^j s^{2k}\frac {\varGamma (k-1/2)}{\varGamma (k+1)}\frac {\varGamma (3/2+j-k)}{2(j-k)+1}\right ] \end{align}

for which the leading, zeroth-order term is

(3.2) \begin{equation} \delta _0(s,t)=2\sqrt {\pi }e^{-s^2/2} \left [\big(1+s^2\big)I_0\big(s^2/2\big)+s^2 I_1\big(s^2/2\big)\right ]-4s, \end{equation}

and $I_0$ and $I_1$ are modified Bessel functions of the first kind. The subsequent, higher-order terms are polynomial and, to fifth-order in the small parameter $1/\lambda$ , the expansion is

(3.3) \begin{equation} \frac {\delta _{in}(s,t)}{\sqrt {\alpha t} f\Delta p/\mu }=\delta _0(s,t)-\frac {1}{\lambda }\frac {4}{\pi }-\frac {1}{\lambda ^3}\frac {2}{3\pi }\big(1-s^2\big)-\frac {1}{\lambda ^5}\frac {1}{5\pi }\big(3-s^2-s^4/2 \big )+O\big(\lambda ^{-7}\big). \end{equation}

3.2. Outer solution

We now examine behaviour over radial distances comparable to the current rupture radius $R(t)$ , relatively much larger than the diffusive length scale $\sqrt {\alpha t}$ . We rescale radial distance as $\rho =r/R(t)$ , and find the asymptotic expansion of the outer solution for slip as $\lambda \rightarrow \infty$ is given by $\delta \sim \delta _{\textit{out}}$ , where

(3.4) \begin{align} \frac {\delta _{\textit{out}}(\rho ,t)}{R(t) f \Delta p/\mu }&=\frac {4}{\pi ^{2}}\sum _{n=0}^\infty \frac {1}{\lambda ^{2n+2}}\frac {1}{\rho ^{2n+1}}\frac {\varGamma (n+3/2)}{2n+1}\frac {\varGamma (n+1/2)}{\varGamma (n+2)} \nonumber\\[3 pt] &\quad\boldsymbol{\cdot }\left [\!\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )+2\sqrt {\pi }\frac {\left(1-\rho ^2\right)^{3/2}}{\rho }\sum _{k=1}^{n}\frac {\rho ^{2k}}{2k+2} \frac {\varGamma (k+2)}{\varGamma (k+1/2)}\right ]\!, \end{align}

where we have deferred discussion of the solution details to Appendix B (and an alternative derivation using a multipole expansion to Appendix C). This expansion, to sixth-order in the small parameter $1/\lambda$ when scaled as such, reduces to

(3.5) \begin{align} \frac {\delta _{\textit{out}}(\rho ,t)}{R(t) f \Delta p/\mu }&=\frac {1}{\lambda ^2}\frac {2}{\pi }\left (\frac {\textrm{acos}(\rho )}{\rho }-\sqrt {1-\rho ^2}\right )+\frac {1}{\lambda ^4}\frac {1}{4 \pi }\left [\frac {\textrm{acos}(\rho )}{\rho ^3}-\sqrt {1-\rho ^2}\left (2-\frac {1}{\rho ^2}\right )\right ] \nonumber \\ &\quad +\frac {1}{\lambda ^6}\frac {3}{16\pi }\left [\frac {\textrm{acos}(\rho )}{\rho ^5}-\sqrt {1-\rho ^2}\left (\frac {8}{3}-\frac {2}{3}\frac {1}{\rho ^2}-\frac {1}{\rho ^4}\right )\right ]+O(\lambda ^{-8}). \end{align}

Figure 2.

Plots of inner (red dashed) and outer (blue dashed) solutions in comparison with full solution (solid black), for a case modestly within the critically stressed regime, $\lambda =2$ . The inner and outer solutions are carried out to (a) $O(\lambda ^{-1})$ , (b) $O(\lambda ^{-3})$ , (c) $O(\lambda ^{-5})$ , (d) $O(\lambda ^{-7})$ and (e) $O(\lambda ^{-9})$ . Panel (e) shows the optimally truncated inner and outer solutions of the expansion (discussion further in § 6.1).

We may compare the inner and outer solutions with direct evaluation of the inversion for slip (2.2). For a fixed value of $\lambda =2$ – only slightly within the critically stressed regime – the panels of figure 2 show the inner and outer solutions as red- and blue-dashed lines, respectively, overlaid on a direct solution for the self-similar slip distribution $\delta$ (solid black line). The series of panels shows the inner and outer solutions truncated to progressively greater order: i.e. the inner and outer series (3.1) and (3.4), when each are scaled by $\sqrt {\alpha t} f\Delta p/\mu$ , are truncated to $O(\lambda ^{-(2N+1)}$ ) with $N$ increasing from 0 to 4 as panels progress from left to right and in descending rows. To within the range examined, the increasing order of the inner and outer asymptotic solutions show progressively better approximation, with the inner solution appearing to capture behaviour over the entire crack. Plotting the ordinate axes of the panels of figure 2 on a logarithmic scale (figure 3), we improve the visible separation of coverage of the inner and outer solutions of their respective regions, as well as the narrowing region of overlap between these solutions.

Figure 3.

Same as in figure 2, except with the ordinate axes now on a logarithmic scale to highlight separation of inner and outer solutions, as well as their diminishing overlap as order is increased to the order of optimal truncation (panel e).

3.3. Matched asymptotics and composite solution

The inner and outer expansions can be used to construct a composite asymptotic solution, which will provide uniform convergence over the entire domain (e.g. Hinch Reference Hinch1991). The composite solution for the critically stressed limit is given by

(3.6) \begin{equation} \frac {\delta _{cs}(\rho ,t)}{\sqrt {\alpha t} f\Delta p/\mu }=\frac {\delta _{in}(\lambda \rho ,t)+\delta _{\textit{out}}(\rho ,t)-\delta _{\textit{overlap}}(\rho ,t)}{\sqrt {\alpha t} f\Delta p/\mu }, \end{equation}

where $\delta _{\textit{overlap}}$ is the common asymptotic behaviour of the inner and outer solutions in an intermediate region corresponding to the outer region ( $s\rightarrow \infty$ ) of the inner solution, and the inner region ( $\rho \rightarrow 0$ ) of the outer solution. This overlap is subtracted to avoid double-counting contributions in this intermediate region when summing the inner and outer solutions.

The expression for the inner solution is given as before, (3.1)–(3.2). Regarding the outer solution, we must carefully note that the scaling for the composite solution (3.6) does not match that of the prior expression of the outer solution (3.4): their scalings differ by a factor of $R(t)/\sqrt {\alpha t}=\lambda$ . Rescaling the outer solution in the manner of (3.6),

(3.7) \begin{align} \begin{split} \frac {\delta _{\textit{out}}(\rho ,t)}{\sqrt {\alpha t} f \Delta p/\mu }&=\frac {4}{\pi ^{2}}\sum _{n=0}^\infty \frac {\varGamma (n+3/2)}{2n+1}\frac {\varGamma (n+1/2)}{\varGamma (n+2)}\frac {1}{(\lambda \rho )^{2n+1}} \\[3 pt] &\quad\boldsymbol{\cdot }\left [\!\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )+2\sqrt {\pi }\frac {\big(1-\rho ^2 \big)^{3/2}}{\rho }\sum _{k=1}^{n}\frac {\rho ^{2k}}{2k+2} \frac {\varGamma (k+2)}{\varGamma (k+1/2)}\right ]\!. \end{split} \end{align}

We may determine the overlap by examining the outer ( $s\rightarrow \infty$ ) behaviour of the inner solution or the inner ( $\rho \rightarrow 0$ ) behaviour of the outer solution. Focusing on the inner solution, its far-field behaviour is given by the power-law decay of the leading order term $\delta _0(s,t)$ , which we determine in Appendix E, plus the polynomial terms contained in the double sum within (3.1) that are the higher-order inner solution terms. These are respectively the first and second sums over $n$ in the following expression for the overlap:

(3.8) \begin{align} \begin{split} \frac {\delta _{\textit{overlap}}(\rho ,t)}{\sqrt {\alpha t} f\Delta p/\mu }=&\frac {2}{\pi }\sum _{n=0}^\infty \frac {\varGamma (n+3/2)}{2n+1}\frac {\varGamma (n+1/2)}{\varGamma (n+2)}\frac {1}{(\lambda \rho )^{2n+1}}\\[3pt] &+\frac {4}{\pi ^2}\sum _{n=0}^\infty \frac {1}{\lambda ^{2n+1}}\frac {1}{2n+1}\left [\sum _{k=0}^n (\lambda \rho )^{2k}\frac {\varGamma (k-1/2)}{\varGamma (k+1)}\frac {\varGamma (3/2+n-k)}{2(n-k)+1}\right ]\!. \end{split} \end{align}

The composite solution is hence

(3.9) \begin{align} \begin{split}& \frac {\delta _{cs}(\rho ,t)}{\sqrt {\alpha t} f\Delta p/\mu }=2 \sqrt {\pi } e^{-(\lambda \rho )^2/2} \left \{[1+(\lambda \rho )^2]I_0\left [\frac {(\lambda \rho )^2}{2}\right ]+(\lambda \rho )^2 I_1\left [\frac {(\lambda \rho )^2}{2}\right ]\right \} -4 \rho \lambda \\[3pt] &\quad+\frac {4}{\pi ^{2}}\sum _{n=0}^\infty \frac {1}{(\lambda \rho )^{2n+1}}\frac {\varGamma (n+3/2)}{2n+1}\frac {\varGamma (n+1/2)}{\varGamma (n+2)}\\[3pt] &\quad\boldsymbol{\cdot } \left [\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )+2\sqrt {\pi }\frac {\big(1-\rho ^2 \big)^{3/2}}{\rho }\sum _{k=1}^{n}\frac {\rho ^{2k}}{2k+2} \frac {\varGamma (k+2)}{\varGamma (k+1/2)}-\frac {\pi }{2}\right ]\!. \end{split} \end{align}

The last term in brackets may be slightly simplified using the identity $\textrm{acos}(\rho )-\pi /2=-\text{asin}(\rho )$ .

To illustrate the asymptotic matching, we examine the above-mentioned series truncated to three terms. In this case, the inner solution is

(3.10) \begin{align} \frac {\delta _{in}(s,t)}{\sqrt {\alpha t}f \Delta p/\mu }=\,\,&2 \sqrt {\pi } e^{-s^2/2} \big [ \big(1+s^2 \big)I_0 \big(s^2/2 \big)+s^2 I_1\big(s^2/2 \big)\big ] -4s \nonumber \\ &-\frac {4}{\pi \lambda }-\frac {2}{3\pi \lambda ^3}\big (1-s^2\big )-\frac {1}{5\pi \lambda ^5}\left (3-s^2-\frac {s^4}{2}\right ) +O\big(\lambda ^{-7} \big) \end{align}

and the outer solution, rescaled in the manner of (3.6), is

(3.11) \begin{align} \frac {\delta _{\textit{out}}(\rho ,t)}{\sqrt {\alpha t} f \Delta p/\mu }=&\frac {1}{\lambda }\frac {2}{\pi }\left (\frac {\textrm{acos}(\rho )}{\rho }-\sqrt {1-\rho ^2}\right )+\frac {1}{\lambda ^3}\frac {1}{4 \pi }\left [\frac {\textrm{acos}(\rho )}{\rho ^3}-\sqrt {1-\rho ^2}\left (2-\frac {1}{\rho ^2}\right )\right ] \nonumber \\ &+\frac {1}{\lambda ^5}\frac {3}{16\pi }\left [\frac {\textrm{acos}(\rho )}{\rho ^5}-\sqrt {1-\rho ^2}\left (\frac {8}{3}-\frac {2}{3}\frac {1}{\rho ^2}-\frac {1}{\rho ^4}\right )\right ]+O(\lambda ^{-8}). \end{align}

We may now evaluate the so-called overlap of the inner and outer solutions, here both expanded to fifth-order in $\lambda ^{-1}$ when scaled as such.

The overlap of these two solutions is determined by comparing the outer behaviour of the inner solution, i.e. the asymptotic behaviour of $\delta _{in}(s\rightarrow \infty ,t)$ with the inner behaviour of the outer solution, i.e. the asymptotic behaviour of $\delta _{\textit{out}}(\rho \rightarrow 0$ ,t). The latter can be determined to be

(3.12)

where we again note here that the scaling of $\delta _{\textit{out}}$ on the left-hand side is different from the scaling of slip in (3.4), in which $R(t)=\lambda \sqrt {\alpha t}$ was used in place of $\sqrt {\alpha t}$ . In (3.12), we scale slip with $\sqrt {\alpha t}$ so that we may correctly compare this inner behaviour with the outer behaviour of the inner solution, which is

(3.13)

Recognising that $s=\rho \lambda$ , we can identify the common, overlapping terms for the two expansions, which we colour in red, blue and green to facilitate comparisons. This overlap is clearly

(3.14)

The composite solution to fifth order in the small parameter $1/\lambda$ is then

(3.15) \begin{align} \frac {\delta _{cs}(\rho ,t)}{\sqrt {\alpha t} f\Delta p/\mu }&=2 \sqrt {\pi } e^{-(\lambda \rho )^2/2} \left \{[1+(\lambda \rho )^2]I_0\left [\frac {(\lambda \rho )^2}{2}\right ]+(\lambda \rho )^2 I_1\left [\frac {(\lambda \rho )^2}{2}\right ]\right \} -4\lambda \rho \nonumber\\[3pt] &\quad-\frac {1}{\lambda }\frac {2}{\pi }\left (\frac {\text{asin}(\rho )}{\rho }+\sqrt {1-\rho ^2}\right )-\frac {1}{\lambda ^3}\frac {1}{4 \pi }\left [\frac {\text{asin}(\rho )}{\rho ^3}+\sqrt {1-\rho ^2}\left (2-\frac {1}{\rho ^2}\right )\!\right ] \nonumber\\[3pt] &\quad-\frac {1}{\lambda ^5}\frac {3}{16\pi }\left [\frac {\text{asin}(\rho )}{\rho ^5}+\sqrt {1-\rho ^2}\left (\frac {8}{3}-\frac {2}{3}\frac {1}{\rho ^2}-\frac {1}{\rho ^4}\right )\right ]+O(\lambda ^{-7}). \end{align}

We now examine the accuracy of the composite solution to represent slip that is only modestly in the critically stressed regime, considering an example case for which $\lambda$ is only moderately large: $\lambda =2$ . In figure 4(a), we compare the composite solution constructed up to ninth order in the small parameter $\lambda ^{-1}$ to a more accurate numerical solution for the complete slip profile. Beyond fifth order, the composite solution is indistinguishable from the numerical solution at this scale. To highlight that the composite solutions do deviate from the precise and accurate numerical solution at various orders, in figures 4(c) and 4(d), we show the discrepancy from the numerical solution on a logarithmic scale at the crack centre $r=0$ and crack tip $r=R(t)$ , respectively. In figure 4(b) (solid lines), we show the distribution along the crack of the absolute error of the composite solution at various orders. We compare this error with the value of the next-order term in the composite solution (dashed lines). We see that by the seventh-order solution (solid cyan), much of the error is contained in the next-order term (underlying dashed cyan). At ninth order, we see a dramatic drop in error, with the next-order term (dotted orange) having a much larger value. This is an indication that, for $\lambda =2$ , the composite solution is most accurate at ninth order. Carrying out the composite solution to higher order will only increase the error. We postpone further discussion of optimal truncation to § 6.1. We see that even for mildly critically stressed conditions, the composite solution is capable of achieving absolute errors of the order of $10^{-5}$ , small in comparison with the $O(1)$ value of slip along the crack.

Figure 4.

(a) Comparison of composite solutions for slip at various orders $\delta _{cs}$ (dashed lines) against the solution for slip $\delta$ (black line) at a value of $\lambda =2$ only modestly in the critically stressed regime. Composite solutions are carried out to the following orders in the small parameter $\lambda ^{-1}$ : leading first (red), third (blue), fifth (green), seventh (cyan) and ninth (orange). By fifth order, the composite solution is indistinguishable from the true solution. (b) Absolute error of composite solutions at various orders are shown as a solid line, whose colours indicate order of composite solution and correspond to those of dashed lines in panel (a). Dashed and dotted lines indicate value of the next-order term in the composite solution: e.g. the red-dashed line shows that the value of the third-order ( $\lambda ^{-3}$ ) term of the composite solution is equal to more than 50 % of the error of the first-order solution (red-solid line). Note the dramatic reduction in error for the composite solution at ninth order (solid orange line), reflecting optimal truncation of the composite solution (discussed further in § 6.1). (c) Deviation of the composite solution $\delta _{cs}$ at various orders (dashed lines) from the solution for slip $\delta$ near the crack centre, $r=0$ . Excellent agreement of optimally truncated composite solution at ninth order is again apparent. (d) Deviation of the composite solution near the crack tip, $r=R(t)$ .

4. Fault slip under marginally pressurised conditions

In the limit of large $T$ (small $\lambda$ ), the self-similar fluid pressure profile stretches faster than crack growth: rupture lags the injected fluid. In this case, we may perform a regular perturbation expansion for the slip profile in powers of the small parameter $\lambda$ . The details of the solution may be found in Appendix D, but the perturbation to arbitrary order is given by

(4.1) \begin{align} \frac {\delta _{\textit{mp}}(\rho ,t)}{R(t)f \Delta p/\mu }&=\frac {8}{\pi }\left (\sqrt {1-\rho ^2}-\rho \text{ acos}(\rho )\right )+2\big (1-\rho ^2\big )^{3/2}\sum _{n=1}^\infty \frac {\lambda ^{2n}}{2n+1}\frac {(-1)^n}{\varGamma (n+3/2)}\nonumber\\&\quad\boldsymbol{\cdot }\sum _{k=1}^{n}\frac {\varGamma (k)}{\varGamma (k+1/2)}\rho ^{2k-2}. \end{align}

To $O(\lambda ^4)$ , the asymptotic solution is

(4.2) \begin{align} \frac {\delta _{\textit{mp}}(\rho ,t)}{R(t)f \Delta p/\mu }&=\frac {8}{\pi }\left (\sqrt {1-\rho ^2}-\rho \text{ acos}(\rho )\right )-\lambda ^2\frac {16}{9\pi } \big(1-\rho ^2 \big)^{3/2}\nonumber\\&\quad+\lambda ^4\frac {32}{225\pi } \big(1-\rho ^2 \big)^{3/2} \big(3+2\rho ^2 \big)+O \big(\lambda ^6 \big). \end{align}

In figure 5, we examine a mildly marginally pressurised solution ( $\lambda =1/2$ ). Panel (a) compares the full solution (solid black curve) with the marginally pressurised solution truncated to various orders: $O(\lambda ^0)$ , $O(\lambda ^2)$ , $O(\lambda ^4)$ and $O(\lambda ^6)$ . Panel (b) shows relative error at each order. In contrast with the critically stressed case, in which the absolute error showed uniform convergence, here, it is the relative error that is uniform, such that the divergence of the type seen in, for example, figure 4(d) does not occur for the truncated marginally pressurised solutions. The regularly decreasing error with increase in order in the series expansions exhibits the standard convergence expected for a regular perturbation expansion. We discuss this convergence further in § 6.2.

Figure 5.

(a) Comparison of the marginally pressurised series solution truncated at various orders (dashed lines) with the full solution (black solid line), both for $\lambda =1/2$ . The series are carried up to: $O(\lambda ^0)$ (red), $O(\lambda ^2)$ (blue), $O(\lambda ^4)$ (green), $O(\lambda ^6)$ (cyan) and $O(\lambda ^8)$ (orange). Beyond the leading order, the marginally pressurised solutions are indistinguishable at this scale. (b) Plot of the relative error of the truncated series solutions, with the same colour scale as panel (a). The regular perturbation series exhibits uniform convergence in terms of relative error.

5. Aseismic moment release; peak slip at injection point

We now calculate two quantities, one characterising far-field behaviour and the other local, as they each depend on $\lambda$ in the marginally pressurised and critically stressed limits. The first quantity is the so-called aseismic moment, which is the leading-order moment in a multipole expansion of the far-field behaviour (e.g. displacement) due to fault slip. However, we also calculate a near-field quantity: the peak slip, which occurs at the crack centre.

The aseismic moment is given by (e.g. Rice Reference Rice1980)

(5.1) \begin{equation} M(t)=\mu \int _0^{R(t)} \delta (r,t) 2\pi r\,{\rm d}r. \end{equation}

For a marginally pressurised fault, we substitute $\delta =\delta _{\textit{mp}}$ in this expression for $M$ and find that

(5.2) \begin{align} \frac {M(t)}{R(t)^3 f\Delta p}&=\frac {16}{9}+\frac {3\pi ^{3/2}}{2}\sum _{n=1}^\infty \frac {\lambda ^{2n}}{2n+1} \frac {(-1)^n}{\varGamma (n+3/2)}\sum _{k=1}^n \frac {\varGamma (k)^2}{\varGamma (k+1/2)\varGamma (k+5/2)}\nonumber \\ &=\frac {16}{9}\left (1- \frac {2}{5}\lambda ^2+\frac {4}{35}\lambda ^4 - \frac {9512}{385\,875}\lambda ^6+ O \big(\lambda ^{8} \big)\right). \end{align}

For a critically stressed fault, we substitue $\delta =\delta _{cs}$ and find that

(5.3) \begin{align} &\frac {M(t)}{R(t)^2 \sqrt {\alpha t} f\Delta p }\sim e^{-\lambda ^2/2}\frac {2\pi ^{3/2}}{3}\left [(3+2\lambda ^2)I_0(\lambda ^2/2)+(1+2\lambda ^2)I_1(\lambda ^2/2)\right ]-\frac {8\pi }{3}\lambda \notag \\&\quad+\frac {8}{\pi }\sum _{n=0}^\infty \frac {1}{\lambda ^{2n+1}}\frac {\varGamma (n+3/2)}{2n+1}\frac {\varGamma (n+1/2)}{\varGamma (n+2)} \boldsymbol{\cdot }\int _0^1 \frac {1}{\rho ^{2n+1}}\left [{-}\left (\text{asin}(\rho )+\rho \sqrt {1-\rho ^2}\right )\right . \notag \\&\quad\left .+\,2\sqrt {\pi }\frac {\left(1-\rho ^2\right)^{3/2}}{\rho }\sum _{k=1}^{n}\frac {\rho ^{2k}}{2k+2} \frac {\varGamma (k+2)}{\varGamma (k+1/2)}\right ] {\rm d}\rho \notag \\&\quad\sim \frac {8}{3}\lambda ^{-1}-\frac {4}{3}\lambda ^{-3}-\frac {2}{3}\lambda ^{-5}-\lambda ^{-7}-\frac {5}{2}\lambda ^{-9} -\frac {35}{4}\lambda ^{-11} +O(\lambda ^{-13}) .\end{align}

In figure 6(a), we plot the variation of moment release with $\lambda$ and compare with the above-mentioned perturbation expansions. In figure 6(a), we normalise the moment release in the manner of (5.2); multiplying (5.3) by a factor of $\lambda ^{-1}$ normalises the moment release in the critically stressed limit in the manner of (5.2). If we instead multiply (5.3) by a factor of $\lambda$ , we find that to leading order, $M(t)\sim (8/3) \lambda (\alpha t)^{3/2} f \Delta p$ .

Figure 6.

Comparison of (a) aseismic moment and (b) peak slip as a function of $\lambda$ (solid black) compared with truncated marginally pressurised (blue dashed) and critically stressed (red dashed) expansions. (a) The marginally pressurised expansion is $O(\lambda ^{8})$ and the critically stressed expansion is $O(\lambda ^{-6})$ . The former reflects the first four terms of (5.2); the latter reflects the first three terms of (5.3) and a $\lambda ^{-1}$ scaling factor due to a change in normalisation. (b) The marginally pressurised expansion $\delta _{\textit{mp}}$ is $O(\lambda ^9)$ and the critically stressed expansion $\delta _{cs}$ is $O(\lambda ^{-5})$ . The former reflects the first five terms of (5.5) and a $\lambda$ scaling factor due to a change in normalisation. The latter reflects the first three terms of (5.4). In both panels, truncation of the critically stressed expansions at three terms is optimal for values of $\lambda$ close to 1. Optimal truncation is discussed in § 6.1.

We may likewise provide asymptotic expressions for the maximum slip of the rupture, $\delta (0,t)$ . For critically stressed conditions, we may obtain the peak slip directly from the inner solution by evaluating (3.1) for $\delta _{in}(0,t)$ ,

(5.4) \begin{align} \frac {\delta (0,t)}{\sqrt {\alpha t} f \Delta p/\mu }&\sim 2\sqrt {\pi }\left (1-\frac {4}{\pi ^2}\sum _{j=0}^\infty \frac {1}{\lambda ^{2j+1}}\frac {\varGamma (3/2+j)}{(1+2j)^2} \right )\nonumber \\ &\sim 2\sqrt {\pi } - \frac {4}{\pi }\frac {1}{\lambda }-\frac {2}{3\pi }\frac {1}{\lambda ^3}-\frac {3}{5\pi }\frac {1}{\lambda ^5}-\frac {15}{14\pi }\frac {1}{\lambda ^7}-\frac {35}{12\pi }\frac {1}{\lambda ^{9}}-\frac {945}{88\pi }\frac {1}{\lambda ^{11}}\nonumber \\ &+O(\lambda ^{-13}). \end{align}

For marginally pressurised conditions, we obtain an expression for the peak slip directly from (4.1),

(5.5) \begin{align} \frac {\delta (0,t)}{R(t) f \Delta p/\mu }&=\frac {4}{\sqrt {\pi }}\sum _{n=0}^\infty \frac {\lambda ^{2n}}{2n+1}\frac {(-1)^n}{\varGamma (n+3/2)} \nonumber \\ &= \frac {8}{\pi }\left (1-\frac {2}{9}\lambda ^2+\frac {4}{75}\lambda ^4-\frac {8}{735}\lambda ^6 + O \big(\lambda ^{8} \big)\right ) .\end{align}

In figure 6(b), we plot the variation of the peak slip with $\lambda$ and compare with the above-mentioned perturbation expansions.

6. Perturbation series error estimation

In the preceding section, carrying out the asymptotic expansion for moment release and peak slip to higher orders makes apparent the different nature of the perturbation expansions under critically stressed and marginally pressurised conditions. For the latter, the expansion is regular and the coefficients in the series diminish in such a manner that the series is convergent, as may be shown by a ratio test. For the former, the perturbation is singular and the series is divergent.

6.1. Optimal truncation of divergent series under critically stressed conditions

In spite of the growing coefficients, we may find an optimal truncation of the divergent series. Specifically, we make use of the uniform convergence, with respect to radial distance $r$ , of the composite solution to freely choose one point along crack to examine the diminishing benefit of keeping additional terms in the series. For convenience, we examine the asymptotic behaviour of slip at the crack centre. We may abstract the expression (5.4) for the slip there as

(6.1) \begin{equation} \frac {\delta (0,t)}{\sqrt {\alpha t} f \Delta p/\mu }\sim 2\sqrt {\pi } + \sum _{j=0}^\infty a_j\lambda ^{-(2j+1)} \end{equation}

and note that the series coefficients diverge for large $j$ in the manner

(6.2) \begin{equation} \frac {a_{j+1}}{a_j}\sim j-\frac {1}{2} +O\big(j^{-1}\big) .\end{equation}

Comparing sequential terms, we are interested in keeping terms provided that, for a given $\lambda$ , the ratio of sequential terms is diminishing, or

(6.3) \begin{equation} \frac {a_{j+1}\lambda ^{-(2j+3)}}{a_j\lambda ^{-(2j+1)}}\lt 1. \end{equation}

Using the leading asymptotic behaviour of the coefficient ratio and rearranging, this is equivalent to the requirement that

(6.4) \begin{equation} j\lt \lambda ^2+1/2. \end{equation}

In other words, to provide the optimal approximation for the peak slip using the series expansion in (5.4) with summation index $j$ (or the optimal approximation of the entire slip profile using the composite solution (3.9) with summation index $n$ ), the series should be carried out to the term whose summation index is the largest integer less than $\lambda ^2$ + 1/2.

We may approximate the error of this optimal truncation in the critically stressed regime. We may approximate this error as the remainder left if we truncate the series (6.1) at a large value of $j$ . Asymptotically, we can approximate this remainder as $a_j$ , with the choice of optimal truncation $j\sim \lambda ^2$ . Recalling that $a_j$ stands in for terms of the series (5.4), this error estimate $a_{\lambda ^2}$ evaluates as

(6.5) \begin{equation} 2\sqrt {\pi }\frac {4}{\pi ^2}\frac {1}{\lambda ^{2 \lambda ^2+1}}\frac {\varGamma \left(3/2+\lambda ^2\right)}{\left(1+2 \lambda ^2\right)^2} ,\end{equation}

which we may simplify using the Stirling-like formula $\varGamma (3/2+x)\sim \sqrt {2\pi } x (x/e)^x$ for large $x$ to

(6.6) \begin{equation} \frac {2^{3/2}}{\pi }\frac {e^{-\lambda ^2}}{\lambda ^3}, \end{equation}

which shows the optimal truncation to be exponentially accurate for large $\lambda$ .

Figure 7.

(a) Plot of absolute error at the crack centre of the asymptotic series solution (5.4) truncated at $j=N$ , for $\lambda =2$ . The optimal truncation is expected to be $N=\lambda ^2$ , and we observe a dramatic drop in error for $N=4$ . (b) Absolute error of the composite solution at the crack centre optimally truncated at $N=\lambda ^2$ , plotted for fixed values of $\lambda$ (black dots). The expected dependence of the error on $\lambda$ (6.6) is shown as a solid line. (Bottom) Spatial distribution of absolute error showing uniform convergence for five composite solutions (3.9) optimally truncated at $N=\lambda ^2$ : $\lambda =1$ (red), $\lambda =2$ (blue), $\lambda =4$ (green), $\lambda =6$ (cyan) and $\lambda =8$ (orange).

We confirm exponential accuracy of the asymptotic approximations in the critically stressed limit by comparisons with analytical or sufficiently accurate numerical solutions. Examining slip at the crack centre, $r=0$ , an analytical expression for slip follows from substituting $s=0$ into (A9),

(6.7) \begin{equation} \frac {\delta (0,t)}{\sqrt {\alpha t} f \Delta p/\mu }=\frac {8}{\pi } \int _0^\lambda \int _0^1 \frac {x e^{-(x\xi )^2}}{\sqrt {1-x^2}}\,{\rm d}x\,{\rm d}\xi =\frac {8}{\pi } \,{}_2 F_2\!\left (\begin{matrix}1/2& 1\\3/2&3/2\end{matrix};-\lambda ^2\right )\lambda, \end{equation}

where ${}_pF_{q}$ is the generalised hypergeometric function. We use the difference between the composite solution and this analytical solution to determine the composite solution’s absolute error. We first evaluate the absolute error of the composite solution at the crack centre for the particular case $\lambda =2$ , truncating the series (3.1) for $\delta _{in}$ at $j=N$ , where $N$ is varied. The results are shown in figure 7(a). The error decreases monotonically until $N$ approaches $\lambda ^2=4$ , at which point, a substantial decrease in error occurs and including further terms in the asymptotic expansion results in a worsening approximation. This optimal series being that truncated with $N=\lambda ^2$ is consistent with the anticipated condition (6.4). We now look to confirm the exponential accuracy of the asymptotic approximations as $\lambda$ increases. To do so, we evaluate the absolute error of the optimally truncated composite solution at the crack centre, truncating the series expression (3.1) for $\delta _{in}$ when the summation index $j=N$ , where $N=\lambda ^2$ . We dub this optimally truncated expansion $\delta _{in,opt}(0,t)$ . We take the difference between this and (6.7) for $\delta (0,t)$ to find the absolute error at fixed values of $\lambda$ , with results shown in figure 7(b) (black dots). Comparing this to the approximated error (6.6) (figure 7 b; red line), we find excellent agreement. Now, to examine the error over the entire crack at fixed values of $\lambda$ , we first truncate the series of the composite solution (3.9), expressing it up to $n=\lambda ^2$ and dubbing this optimal solution $\delta _{cs,opt}$ . Subtracting this from the slip $\delta$ , we obtain the absolute error across the crack. Repeatedly calculating this error at fixed values of $\lambda$ , we find that the error first calculated at the crack centre (figure 7 b) is essentially uniform across the crack (figure 7 c), demonstrating uniform convergence and exponential accuracy of the composite solution.

Figure 8.

(a) Absolute error of the slip at the crack centre for the marginally pressurised series solution at $O(\lambda ^{2n})$ for the particular case $\lambda =1/2$ (black dots). The expected error (6.10) is shown as a solid line. (b) Distribution of the absolute error along the expanding crack for the case $\lambda =1/2$ truncated at $O(\lambda ^0)$ (red), $O(\lambda ^2)$ (blue), $O(\lambda ^4)$ (green), $O(\lambda ^6)$ (cyan), $O(\lambda ^8)$ (orange).

6.2. Convergence of regular perturbation under marginally pressurised conditions

To estimate the absolute error of truncating the series expansion of the solution under marginally pressurised conditions, we begin with the expression for the peak slip (5.5), abstracted to

(6.8) \begin{equation} \frac {\delta (0,t)}{R(t) f \Delta p/\mu }= \sum _{j=0}^\infty b_j\lambda ^{2j}. \end{equation}

Given that the series is a regular, convergent expansion, we may estimate the error of truncating this series at $j=n$ , by the value of the $j=(n+1)$ th term. We may approximate this term for large $n$ , noting that in this limit,

(6.9) \begin{equation} b_{n+1}\sim \frac {\sqrt {2}}{\pi }\frac {e^n}{n^{n+3}}, \end{equation}

such that the error estimate is

(6.10) \begin{equation} \frac {\sqrt {2}}{\pi }\frac {e^n}{n^{n+3}}\lambda ^{2n+2}. \end{equation}

In figure 8(a), we calculate the absolute error incurred by truncating the series (5.5) at $j=n$ by subtracting this value from the analytical expression (6.7), and comparing this error with the above-mentioned estimate. We confirm that, even for conditions only modestly in the marginally pressurised regime, the value of the error decreases rapidly with an increasing number of terms. Indeed, keeping only two terms in the expansion is sufficient for an absolute error of $10^{-3}$ . In figure 8(b), we plot the absolute error over the extent of the crack for various truncations of the series of the marginally pressurised series solution, $\delta _{\textit{mp}}$ , where we see the absolute error is greatest at the crack centre. From figure 5(b), in which we previously plotted the relative error for the same mildly pressurised case ( $\lambda =1/2$ ), we found that the relative error is greatest at the crack tip, has a maximum of approximately 18 % for the zeroth-order solution, and dramatically decreases to 1 % or 0.1 % when keeping one or two additional terms (i.e. up to second- or fourth-order in $\lambda$ ), respectively.

7. Discussion

7.1. Choice of perturbation parameter

In the preceding, the perturbation parameter has been $\lambda$ , the ratio of the crack radius $R$ to the diffusion length scale $\sqrt {\alpha t}$ . This may be a quantity one may approximately infer from field measurements, say by seismological observation of seismicity migration suspected to be caused by fluid-induced aseismic slip, in which case, one could in turn infer the corresponding value of initial state of fault-resolved tractions and strength of injection, reflected in the parameter $T$ . However, there may be instances when $\lambda$ is not a quantity one would necessarily expect to know a priori. Rather, if one is given field measurements or inferences of stress state and fluid pressure, one would instead be in a position to calculate a value for $T$ and then look for the implied value of $\lambda$ . How might one do so? A straightforward route would be to estimate $\lambda$ , given $T$ , by tabulating numerical evaluations of the asymptotic expansions (2.19) or (2.21), which can be truncated to sufficiently match the true solution over the entire parameter range (e.g. as in figure 1). This would circumvent a direct numerical solution of (2.13) or a cumbersome inversion of the hypergeometric series of Sáez et al. (Reference Sáez, Lecampion, Bhattacharya and Viesca2022). However, when looking for expressions of slip and moment release, $\lambda$ is the natural perturbation parameter: given the initial problem formulation, one may directly proceed to find the asymptotic behaviour for $T$ and slip $\delta$ as series expansions in powers of $\lambda$ , as demonstrated.

One may nonetheless ask whether asymptotic expansions could be conducted in terms of $T$ rather than $\lambda$ . Of the two quantities, $\lambda$ and $\delta$ , it is only necessary to look for such expansions of the former: once a value of $\lambda$ is approximated for a given $T$ , one may proceed with the provided expansions for slip and moment release in terms of $\lambda$ to the desired accuracy. With this in mind, we provide asymptotic expansions of $T$ in terms of $\lambda$ in the following.

In the critically stressed limit, we may look for an asymptotic expansion of $\lambda$ as the power series $\lambda \sim \sum _{n=-1}^\infty A_{2n+1} u^{2n+1}$ , where $u=\sqrt {2T}$ , which we may in turn substitute for $\lambda$ in (2.13),

(7.1) \begin{equation} T=\int _0^1 \frac {x\, E_1\!\left [\left (\displaystyle {A_{-1}}{u^{-1}}+A_1 u +A_3 u^3+A_5 u^5 \ldots \right )^2 x^2\right ]}{\sqrt {1-x^2}}\,{\rm d}x. \end{equation}

After asymptotically approximating the integral in the limit $u\rightarrow 0$ (i.e. $T\rightarrow 0$ ), we may successively solve for the coefficients $A_n$ by matching each side of (7.1) at increasing order of $u$ . To the leading four orders, this reduces to

(7.2) \begin{align} O(u^2)\!: \quad 1&=\frac {1}{A_{-1}}; \nonumber \\ O(u^4)\!: \quad 0&=1-8A_{-1}A_1;\nonumber \\ O(u^6)\!: \quad 0&=1-4 \big (A_{-1}A_1-3A_{-1}A_{1}^2+2A_{-1}^3 A_3\big );\nonumber \\ O(u^8)\!: \quad 0&=15-16\big (3A_{-1}A_1-5A_{-1}^2A_1^2+8 A_{-1}^3A_1^3 \nonumber \\ &\quad +\,2A_{-1}^3A^{}_3-12A_{-1}^4A^{}_1A^{}_3+4A_{-1}^5A^{}_5\big ), \end{align}

from which we find that $A_{-1}=1$ , $A_1=1/8$ , $A_3=11/128$ and $A_5=149/1024$ providing for the expansion of $\lambda$ in the critically stressed limit, $T\rightarrow 0$ ,

(7.3) \begin{equation} \lambda \sim \frac {1}{\sqrt {2T}}+\frac {1}{8}\sqrt {2T} +\frac {11}{128}\sqrt {2T}^3+\frac {149}{1024}\sqrt {2T}^5+O \big(T^{7/2} \big) \end{equation}

a divergent series, like the expansion of $T$ in terms of $\lambda$ in the limit of large $\lambda$ .

In the marginally pressurised limit, we may look for an expansion of $\lambda$ as the series $\lambda =\sum _{n=0}^\infty B_{2n+1}w^{2n+1}$ , where $w=(1/2)\exp [(2-\gamma -T)/2]$ , to be substituted for $\lambda$ in (2.13),

(7.4) \begin{align} T=\int _0^1 \frac {x\, E_1\!\left [\left (B_1 w+B_3w^3+B_5 w^5 +\cdots \right )^2 x^2\right ]}{\sqrt {1-x^2}}\,{\rm d}x. \end{align}

This time, we may approximate the integral in the limit $w\rightarrow 0$ (i.e. $T\rightarrow \infty$ ) and successively solve for the coefficients $B_n$ by matching each side of (7.4) at increasing order of $w$ . To the leading four orders, this reduces to

(7.5) \begin{align} O(\log w)\!: & \quad 0=\log (B_1); \nonumber \\ O(w^2)\!:& \quad 0=B_1^3-3B_3 ;\nonumber \\ O(w^4)\!:& \quad 0=2B_1^6-5\big (4B_1^3B_3+3 B_3^2-6B_1B_5\big );\nonumber \\ O(w^6)\!:& \quad 0=4B_1^9-84B_1^6B_3+105 \big (B_1^3B_3^2-B^3_3+2B_1^4B_5+3B_1B_3B_5-3B^2_1B_7\big ), \end{align}

from which we find that $B_1=1$ , $B_3=1/3$ , $B_5=19/90$ and $B_7=181/1134$ , providing for the regular expansion of $\lambda$ in the marginally pressurised limit, $T\rightarrow \infty$ ,

(7.6) \begin{equation} \lambda = \frac {1}{2}e^{\frac {1}{2}(2-\gamma -T)}+\frac {1}{24}e^{\frac {3}{2}(2-\gamma -T)}+\frac {19}{2880}e^{\frac {5}{2}(2-\gamma -T)}+\frac {181}{145152}e^{\frac {7}{2}(2-\gamma -T)}+O\left(e^{-\frac {9}{2}T}\right) \end{equation}

Figure 9.

Plot of the variation of the prefactor $\lambda$ for the self-similar expansion of the crack radius, $R(t)=\lambda \sqrt {\alpha t}$ (solid black), as it depends on the stress-injection parameter $T$ , which characterises the pre-injection state of fault stress and strength as well as the strength of the injection. This plot is the inverse of figure 1 and showcases the asymptotic expansions of $\lambda$ in terms of $T$ : in comparison are the critically stressed (7.3, red dashed) and the marginally pressurised (7.6, blue dashed) perturbation expansions in the small and large $T$ limit, respectively. The marginally pressurised expansion shown here is carried to $O(e^{(-7/2)T})$ and the critically stressed expansion is carried to $O(T^{(1/2)})$ , which represents an optimal truncation for $T$ near 1.

7.2. Why $\nu =0?$

The particular case studied here, in which Poisson’s ratio $\nu =0$ , may seem a peculiar, unphysical basis; however, this solution can be seen as the more readily solvable leading-order asymptotic solution for the more general problem of $\nu \neq 0$ . Rupture propagation when $\nu \neq 0$ is no longer axisymmetric and a free-boundary problem needs to be resolved, determining the boundary shape of the crack. This point was discussed in the preceding work of Sáez et al. (Reference Sáez, Lecampion, Bhattacharya and Viesca2022), who used a boundary-element method to explore solutions for $\nu$ in the range $0\leqslant \nu \leqslant 0.45$ . They observed in their numerical solutions that the growth of rupture appears to remain self-similar when $\nu \neq 0$ , which itself may be unsurprising considering that the only length scale in the problem is $\sqrt {\alpha t}$ . Moreover, they observed that the self-selected shape of the expanding rupture was well approximated by an ellipse over much of the parameter space, denoting the expanding length of the major and minor axes as $2a(t)$ and $2b(t)$ , respectively. Remarkably, they observed that, for a given value of $T$ , the area of this self-similarly expanding ellipse agreed precisely with the area of the self-similar expanding circular crack that occurs for $\nu =0$ : i.e. to within an error of 1 % or less, they observed that $\pi R(t)^2=\pi a(t) b(t)$ for a given $T$ and $\sqrt {\alpha t}$ , independent of the value of $\nu$ . This permitted the generalisation of the pre-factor $\lambda$ from the $\nu =0$ definition ( $\lambda = R(t)/\sqrt {\alpha t}$ ) to the cases where $\nu \neq 0$ , ( $\lambda = \sqrt {a(t)b(t)}/\sqrt {\alpha t}$ ): the solution of $\lambda (T)$ from the axisymmetric $\nu =0$ was precisely that $\lambda (T)$ for the non-axisymmetric case. In addition, Sáez et al. (Reference Sáez, Lecampion, Bhattacharya and Viesca2022) observed in their numerical solutions that the ratio of major to minor axes behaved asymptotically as $a/b\sim (1-\nu )^{-1}$ and $(3-\nu )/(3-2\nu )$ in the critically stressed and marginally pressurised limits, respectively, over the explored range of $0.05\leqslant \nu \lt 0.45$ .

Thus, Sáez et al. (Reference Sáez, Lecampion, Bhattacharya and Viesca2022) demonstrated that the circular crack solution adequately describes the rupture area and extent, to within a $\nu$ -dependent shape correction in these limits. That and subsequent work has yet to examine, however, how the distribution of slip within the elliptic rupture ( $\nu \neq 0$ ) differs from that of a circular one ( $\nu =0$ ) at the same $T$ , e.g. whether a mapping of the axisymmetric slip from a circular to elliptical domain would suffice. Subsequent work has, however, indicated that the effect of $\nu \neq 0$ on the slip is sufficiently small such that the $\nu =0$ solution can be considered an adequate leading-order solution in a perturbation expansion of the solution about $\nu =0$ . Specifically, Sáez et al. (Reference Sáez, Lecampion, Bhattacharya and Viesca2022) numerically solved for the total moment released following the arrest of the self-similarly expanding aseismic rupture, due to a sudden arrest of an otherwise constant injection rate. Moment release is a useful global measure of the slip as, to within a factor of the shear modulus, it is defined as the integral of slip over the rupture area (see § 5). They found that the moment released for a non-circular rupture ( $\nu =0.25$ ) was consistently 14 % less than circular rupture ( $\nu =0$ ) at the same value of $T$ , explored numerically over two orders of magnitude of $T$ , again implying a $\nu$ -dependent correction to the $\nu =0$ solution.

Thus, the axisymmetric, $\nu =0$ solutions represent a fairly accurate representation of the more general $\nu \neq 0$ problem that would otherwise need to be solved numerically using a boundary-element method for a two-dimensional crack in a three-dimensional medium. Furthermore, the asymptotic solutions presented here provide closed-form expressions of the solutions essentially over the entire range of $T$ , as evidenced, for example, in figures 1 and 6. The utility of an analytical, reduced-order model, for example, is apparent when considering it as an elementary model for the inference of aseismic slip due to a known or suspected, natural or artificial, fluid source that drives aseismic slip. Not having to run an elaborate forward numerical model at each iteration of an inversion drastically reduces computation time. Such an axisymmetric representation was used, for example, by Bhattacharya & Viesca (Reference Bhattacharya and Viesca2019) in the inversion of hydromechanical properties for a time-variable injection history. A constant injection rate history is a reasonable starting point for a case in which data are too poor to constrain an injection history, as would likely be the case for any naturally occurring aseismic rupture in an interplate or crustal fault setting, such as, for example, the observed concentric migration of microseismicity thought to be due to fluid migration (Ross et al. Reference Ross, Cochran, Trugman and Smith2020). Here, for an assumed set of properties of the fault and surrounding rock (friction coefficient, hydraulic diffusivity, shear modulus) and values of the initial fault tractions and quantities related to injection (as they enter lumped sets of parameters determining $\Delta p$ and $T$ ), the aseismic source time function, or moment history $M(t)$ , which enter the far-field representations of quasi-static displacement and stress fields used for comparison with observations, is directly determined by the presented solutions, eliminating the need for a forward-model solution step.

The additional utility of such asymptotic expansions is as a benchmark for numerical methods capable of solving the posed problem and as a basis for further analytical work. The problem explored here is the first instance of a crack in a three-dimensional medium (fluid-driven or otherwise) having a perturbation solution carried to arbitrary order. Sáez et al. (Reference Sáez, Lecampion, Bhattacharya and Viesca2022) were able to use an incomplete leading-order solution to this end, missing only the inner, boundary-layer solution in the critically stressed limit. With the complete perturbation expansions, precise determination of the absolute error of numerical procedures can be attained. This contribution is timely as long-standing initiatives for benchmarking models of sequences of earthquakes and seismic slip (SEAS) (e.g. Harris et al. Reference Harris2009; Erickson et al. Reference Erickson2020) have recently turned their attention towards models of fluid-induced aseismic slip (Lambert et al. Reference Lambert2025). This initiative has incorporated earlier results for simpler, two-dimensional problems (Viesca Reference Viesca2021) in benchmark determination. The methods used here and in previous work (Viesca Reference Viesca2021; Sáez et al. Reference Sáez, Lecampion, Bhattacharya and Viesca2022) have motivated the development of analytical results accounting for more realistic descriptions of fault behaviour, including slip-driven enhancement of the along-fault permeability and fault dilatancy (Dunham Reference Dunham2024). Additionally, recent results also show that such relatively simple models can effectively represent the behaviour, including slip distribution and rupture extent, of models with more elaborate, rate- and state-dependent models for friction coefficient and permeability evolution (Berrios-Rivera, Ozawa & Dunham Reference Berrios-Rivera, Ozawa and Dunham2025).

8. Conclusion

In this article, we provided asymptotic perturbation expansions to arbitrary order for a problem of an expanding, mixed-mode circular shear fracture. We demonstrated that expressing these solutions to arbitrary order allows for accurate representation of the full solution over the entire range of the problem parameter, including when the parameter is neither large nor small. The specific problem pursued here provided the opportunity to illustrate both regular and singular perturbations methods as applied to singular integral equations intrinsic to fracture mechanics, with the singular perturbation involving a composite solution consisting of a boundary layer asymptotically matched with an outer solution. We detailed the derivation of these solutions and, in the case of the outer solution, we provide an alternative derivation by multipole expansion.

Funding

This work was supported by the National Science Foundation (grant number 1653382).

Declaration of interests

The authors report no conflict of interest.

Appendix A. Asymptotic solution details of inner solution for critically stressed conditions

Here, we derive the boundary-layer solution for slip near $r=0$ valid over distances comparable to the diffusion length scale $\sqrt {\alpha t}$ . However, we note that in our original inversion for slip, (2.2) and (2.5), radial distances are scaled by the current rupture radius $R$ . Here, we revisit the problem of solving for slip, but now scale radial distances by $\sqrt {\alpha t}$ , denoting

(A1) \begin{align} s=r/\sqrt {\alpha t}. \end{align}

Substituting the variable of integration in (2.5) with $\eta =\xi /\lambda$ , we may thus rewrite (2.5) as

(A2) \begin{equation} \delta (s,t)=-\frac {4}{\pi }\frac {\sqrt {\alpha t}}{\mu }\int _s^\lambda \sqrt {\xi ^2-s^2}\frac {{\rm d}f}{{\rm d}\xi }\,{\rm d}\xi, \end{equation}

where

(A3) \begin{equation} \frac {{\rm d}f}{{\rm d}\xi }=\frac {d}{{\rm d}\xi }\int _0^1 \frac {x \Delta \tau \big(x\xi \sqrt {\alpha t} \big)}{\sqrt {1-x^2}}\,{\rm d}x, \end{equation}

in which we have suppressed explicit time dependence of $\Delta \tau$ , except through the term $\sqrt {\alpha t}$ . That we may do so is apparent considering that the particular expression for $\Delta \tau$ of interest follows the combination of (2.3) and (2.7)–(2.12) and the requirement that $\tau =\tau _s$ ,

(A4) \begin{align} \Delta \tau (r,t)=\tau _b-f \big [\sigma '-\Delta p E_1 \big(r^2/\alpha t \big)\big ], \end{align}

evaluating

(A5) \begin{align} \Delta \tau (x\xi \sqrt {\alpha t},t)=\tau _b-f\big (\sigma '-\Delta p E_1\big[(x\xi )^2 \big]\big ). \end{align}

In this case, we can simplify $\Delta \tau (x\xi \sqrt {\alpha t},t)\Rightarrow \Delta \tau (x\xi )$ and do so in our evaluation of (A3),

(A6) \begin{align} \frac {{\rm d}f}{{\rm d}\xi }&=\frac {d}{{\rm d}\xi }\int _0^1 \frac {x \Delta \tau (x\xi )}{\sqrt {1-x^2}}\,{\rm d}x\nonumber \\ &=\int _0^1\frac {x}{\sqrt {1-x^2}} \displaystyle \frac {{\rm d}\Delta \tau }{{\rm d}w}{\displaystyle \frac {{\rm d}w}{{\rm d}\xi }}\,{\rm d}x\nonumber \\ &=-2f\Delta p \, \frac {1}{\xi }\int _0^1 \frac {x e^{-(x\xi )^2}}{\sqrt {1-x^2}}\,{\rm d}x, \end{align}

where in (A6), $w=x \xi$ and ${\rm d}\Delta \tau (w)/{\rm d}w=2 \exp [-(x \xi )^2]/(x\xi )$ . We now label the integral

(A7) \begin{equation} F(\xi )=\int _0^1 \frac {x e^{-(x\xi )^2}}{\sqrt {1-x^2}}\,{\rm d}x \end{equation}

such that (A2) becomes

(A8) \begin{align} \delta (s,t)=\frac {\sqrt {\alpha t}}{\mu } f \Delta p\frac {8}{\pi } \int _s^\lambda \sqrt {\xi ^2-s^2}\frac {F(\xi )}{\xi }\,{\rm d}\xi ,\end{align}

or more simply,

(A9) \begin{equation} \frac {\delta (s,t)}{\sqrt {\alpha t} f \Delta p/\mu }=\frac {8}{\pi } \int _s^\lambda \sqrt {1-(s/\xi )^2}\,F(\xi )\,{\rm d}\xi. \end{equation}

We now look to asymptotically evaluate (A9), which is exact, in the limit of large $\lambda$ to arrive to our expansion approximation of the inner solution for critically stressed conditions. To do so, we will first rewrite (A9) as

(A10) \begin{equation} \frac {\delta (s,t)}{\sqrt {\alpha t} f \Delta p/\mu }=\frac {8}{\pi }\left [\int _s^\infty \sqrt {1-(s/\xi )^2}F(\xi )\,{\rm d}\xi -\int _\lambda ^\infty \sqrt {1-(s/\xi )^2}F(\xi )\,{\rm d}\xi \right ]\!, \end{equation}

where the leading-order term, $O(1)$ , is

(A11) \begin{equation} \frac {8}{\pi }\int _s^\infty \sqrt {1-(s/\xi )^2}F(\xi )\,{\rm d}\xi =2\sqrt {\pi }e^{-s^2/2} \big [\big(1+s^2 \big)I_0 \big(s^2/2 \big)+s^2 I_1 \big(s^2/2 \big)\big ]-4s, \end{equation}

where $I_0$ and $I_1$ modified Bessel functions of the first kind. We turn our attention to the second integral in (A10) and note that, for large $\lambda$ , the integration variable $\xi$ takes on only large values warranting a Taylor expansion of the integrand’s square-root term,

(A12) \begin{equation} -\int _\lambda ^\infty \sqrt {1-(s/\xi )^2}F(\xi )\,{\rm d}\xi \approx -\int _\lambda ^\infty \left (1-\frac {1}{2}(s/\xi )^2-\frac {1}{8}(s/\xi )^4+\cdots \right ) F(\xi )\,{\rm d}\xi. \end{equation}

We look to evaluate the series of integrals of (A12), which are the remaining $O(1/\lambda )$ terms of (A10).

We start with the first integral of the expansion (A12),

(A13) \begin{align} I(\lambda )=-\int _\lambda ^\infty F(\xi )\,{\rm d}\xi. \end{align}

To expand $I$ for large $\lambda$ , we note that

(A14) \begin{align} \frac {{\rm d}I}{{\rm d}\lambda }=F(\lambda ) \end{align}

and seek the expansion of

(A15) \begin{equation} F(\lambda )=\int _0^1 \frac {x e^{-(\lambda x)^2}}{\sqrt {1-x^2}}\,{\rm d}x .\end{equation}

This will be useful not only to arrive to $I(\lambda )$ , but also to the remaining integrals in the expansion (A12). We use Laplace’s method for asymptotic expansion, substituting the binomial-like expansion for the Maclaurin series of $(1-x^2)^{-1/2}$ ,

(A16) \begin{equation} \frac {1}{\sqrt {1-x^2}}=\sum _{k=0}^\infty \begin{pmatrix} -1/2\\ k \end{pmatrix} (-x^2)^k=\sum _{k=0}^\infty \frac {(2k)!}{(k!)^2}\left (\frac {x}{2}\right )^{2k}, \end{equation}

and passing the upper limit of integration to $\infty$ without loss of generality for the power-law asymptotic expansion,

(A17) \begin{align} F(\lambda )\sim \sum _{k=0}^\infty \frac {(2k)!}{(k!)^2} \frac {1}{2^{2k}} \int _0^\infty x^{2k+1} e^{-(\lambda x)^2}\,{\rm d}x, \end{align}

we may evaluate the integrals using the change of variable $u=(\lambda x)^2$ ,

(A18) \begin{align} \int _0^\infty x^{2k+1} e^{-(\lambda x)^2}\,{\rm d}x=\frac {1}{2\lambda } \int _0^\infty \left (\frac {\sqrt {u}}{\lambda }\right )^{2k+1} e^{-u}\,{\rm d}u=\frac {\varGamma (k+1)}{2 \lambda ^{2k+2}}=\frac {k!}{2 \lambda ^{2k+2}}, \end{align}

such that

(A19) \begin{align} F(\lambda )\sim \sum _{k=0}^\infty \frac {1}{\lambda ^{2k+2} }\frac {1}{2^{2k+1}} \frac {(2k)!}{k!} .\end{align}

Since $(2k)!=(2k+1)!/(2k+1)$ and $n!=\varGamma (n+1)$ , we may rewrite (A19) as

(A20) \begin{align} F(\lambda )\sim \sum _{k=0}^\infty \frac {1}{\lambda ^{2k+2} }\frac {1}{2^{2k+1}} \frac {\varGamma (2k+2)}{\varGamma (k+1)}\frac {1}{2k+1}. \end{align}

Using Legendre’s duplication formula for the Gamma function

(A21) \begin{equation} 2^{1-m}\,\varGamma (m)/\varGamma (m/2)=\varGamma (m/2+1/2)/\sqrt {\pi } ,\end{equation}

while identifying $m=2k+2$ , we may simplify the expression for $F(\lambda )$ to

(A22) \begin{equation} F(\lambda )\sim \sum _{k=0}^\infty \frac {1}{\sqrt {\pi }} \frac {1}{\lambda ^{2k+2} }\frac {\varGamma (k+3/2)}{2k+1} \end{equation}

and we may finally asymptotically approximate the integral

(A23) \begin{align} I(\lambda )=-\int _\lambda ^\infty F(\lambda )\,{\rm d}\lambda \sim -\frac {1}{\sqrt {\pi }}\sum _{n=0}^\infty \frac {1}{\lambda ^{2n+1}} \frac {\varGamma (n+3/2)}{(2n+1)^2} .\end{align}

We now look to evaluate all terms in the expansion (A12), which can be written in summation form as

(A24) \begin{align} -\int _\lambda ^\infty \sqrt {1-(s/\xi )^2}F(\xi )\,{\rm d}\xi =-\int _\lambda ^\infty \left [-\sum _{k=0}^\infty \frac {(k-3/2)!}{2\sqrt {\pi } k!}\left (\frac {s}{\xi }\right )^{2k}\right ] F(\xi )\,{\rm d}\xi ,\end{align}

where the summation in brackets is the Taylor series of $\sqrt {1-(s/\xi )^2}$ about $s/\xi =0$ . Again passing from factorials to Gamma functions, via the identity $\varGamma (n+1)=n!$ , and rearranging, we have

(A25) \begin{equation} -\int _\lambda ^\infty \sqrt {1-(s/\xi )^2}F(\xi )\,{\rm d}\xi = \sum _{k=0}^\infty \frac {1}{2\sqrt {\pi }}\frac {\varGamma (k-1/2)}{\varGamma (k+1)}s^{2k} \int _\lambda ^\infty \frac {F(\xi )}{\xi ^{2k}} \,{\rm d}\xi, \end{equation}

we denote the integral

(A26) \begin{align} J(\lambda )=\int _\lambda ^\infty \frac {F(\xi )}{\xi ^{2k}}\, {\rm d}\xi \end{align}

such that

(A27) \begin{align} \frac {{\rm d}J}{{\rm d}\lambda }=-\frac {F(\lambda )}{\lambda ^{2k}}\sim -\frac {1}{\sqrt {\pi }}\frac {1}{\lambda ^{2k}}\left (\sum _{m=1,3,\,\ldots }^\infty \frac {1}{\lambda ^{2n+2}}\frac {\varGamma (n+3/2)}{2n+1}\right ) ,\end{align}

in which we have substituted the previously found asymptotic expansion for $F(\lambda )$ , (A22). We may subsequently integrate to find

(A28) \begin{align} J(\lambda )=-\int _\lambda ^\infty \frac {{\rm d}J}{{\rm d}\lambda }\,{\rm d}\lambda \sim \frac {1}{\sqrt {\pi }}\sum _{n=0}^\infty \frac {1}{\lambda ^{2(n+k)+1}}\frac {1}{2(n+k)+1}\frac {\varGamma (n+3/2)}{2n+1}. \end{align}

Substituting this expression for $J$ into (A25),

(A29) \begin{align}& -\int _\lambda ^\infty \sqrt {1-(s/\xi )^2}F(\xi )\,{\rm d}\xi\nonumber\\&\quad\sim \frac {1}{2\pi }\sum _{k=0}^\infty \sum _{n=0}^\infty s^{2k}\frac {\varGamma (k-1/2)}{\varGamma (k+1)}\frac {1}{\lambda ^{2(n+k)+1}}\frac {1}{2(n+k)+1}\frac {\varGamma (n+3/2)}{2n+1} \\[-12pt] \nonumber \end{align}
(A30) \begin{align} &\quad\sim \frac {1}{2\pi }\sum _{k=0}^\infty \sum _{j=k}^\infty s^{2k}\frac {\varGamma (k-1/2)}{\varGamma (k+1)}\frac {1}{\lambda ^{2j+1}}\frac {1}{2j+1}\frac {\varGamma (3/2+j-k)}{2(j-k)+1} \\[-12pt] \nonumber \end{align}
(A31) \begin{align} &\quad\sim \frac {1}{2\pi }\sum _{j=0}^\infty \frac {1}{\lambda ^{2j+1}}\frac {1}{2j+1}\left [\sum _{k=0}^j s^{2k}\frac {\varGamma (k-1/2)}{\varGamma (k+1)}\frac {\varGamma (3/2+j-k)}{2(j-k)+1}\right ]\!, \\[9pt] \nonumber \end{align}

where, in passing from (A29) to (A30), we replaced the summation index $n$ with another $j$ chosen such that $2j+1=2(n+k)+1$ , or $n=j-k$ , and in passing from (A30) to (A31), we swapped the order of summation. We thus find that the remaining term in (A10) is asymptotically expanded for large $\lambda$ as

(A32) \begin{align} \frac {8}{\pi }\left (-\int _\lambda ^\infty \sqrt {1-(s/\xi )^2}F(\xi )\,{\rm d}\xi \right )&\sim -\frac {1}{\lambda }\frac {4}{\pi }-\frac {1}{\lambda ^3}\frac {2}{3\pi }\big(1-s^2 \big)-\frac {1}{\lambda ^5}\frac {1}{5\pi }\big (3-s^2-s^4/2\big )\nonumber\\&\quad+O\big(\lambda ^{-7}\big) .\end{align}

We therefore arrive to the asymptotic expansion of the inner solution for slip in the critically stressed limit, $\delta \sim \delta _{in}$ as $\lambda \rightarrow \infty$ , where

(A33) \begin{equation} \frac {\delta _{in}(s,t)}{\sqrt {\alpha t} f\Delta p/\mu }=\delta _0(s,t)-\frac {1}{\lambda }\frac {4}{\pi }-\frac {1}{\lambda ^3}\frac {2}{3\pi }\big(1-s^2\big)-\frac {1}{\lambda ^5}\frac {1}{5\pi }\big (3-s^2-s^4/2\big )+O \big(\lambda ^{-7} \big) \end{equation}

and

(A34) \begin{align} \delta _0(s,t)=2\sqrt {\pi }e^{-s^2/2} \big [(1+s)I_0\big(s^2/2\big)+s^2 I_1\big(s^2/2\big)\big ]-4s, \end{align}

the outer behaviour of which, as $s\rightarrow \infty$ , is

(A35) \begin{align} \frac {\delta _{in}(s\rightarrow \infty ,t)}{\sqrt {\alpha t} f\Delta p/\mu }&=\left (\frac {1}{s}+\frac {1}{8 s^3}+\frac {3}{32 s^5}+O\big(s^{-7}\big)\right )-\frac {1}{\lambda }\frac {4}{\pi }-\frac {1}{\lambda ^3}\frac {2}{3\pi }\big(1-s^2\big)\nonumber\\&\quad-\frac {1}{\lambda ^5}\frac {1}{5\pi }\big (3-s^2-s^4/2\big )+O \big(\lambda ^{-7} \big) .\end{align}

Appendix B. Asymptotic solution details of outer solution for critically stressed conditions

Here, we turn our attention to the solution for slip for radial distances $r$ comparable to the crack-tip radius $R(t)$ when $R(t)\gg \sqrt {\alpha t}$ , i.e. when $\lambda$ is large. We return to the inversion for slip (2.6), applying the non-singular condition $f(1)=0$ :

(B1) \begin{equation} \delta ({\rho },t)=-\frac {4}{\pi }\frac {R}{\mu }\int _\rho ^1 \sqrt {\eta ^2-\rho ^2}\frac {\rm d}{{\rm d}\eta }\int _0^1 \frac {x \Delta \tau (x\eta R,t)}{\sqrt {1-x^2}}\,{\rm d}x{\rm d}\eta, \end{equation}

where we recall that

(B2) \begin{align} \Delta \tau (r,t)=\tau _b-f \big[\sigma '-\Delta p E_1 \big(r^2/\alpha t \big) \big]. \end{align}

When evaluating

(B3) \begin{align} \Delta \tau (x\eta R(t),t)=\tau _b-f\big (\sigma '-\Delta p E_1 \big[(x\eta \lambda )^2 \big]\big ), \end{align}

we note that we can simplify the left-hand side $\Delta \tau (x\xi \sqrt {\alpha t},t)\Rightarrow \Delta \tau (x\eta \lambda )$ and do so in our evaluation of (B1), as was done similarly for the inner solution in Appendix A:

(B4) \begin{align} \frac {\rm d}{{\rm d}\eta }\int _0^1 \frac {x \Delta \tau (x\eta \lambda )}{\sqrt {1-x^2}}\,{\rm d}x&=\int _0^1 \frac {x}{\sqrt {1-x^2}}\frac {{\rm d}\Delta \tau }{{\rm d}w}\frac {{\rm d}w}{{\rm d}\zeta }\,{\rm d}x\nonumber \\ &=-2f\Delta p \frac {1}{\eta }\int _0^1\frac {xe^{-(x\eta \lambda )^2}}{\sqrt {1-x^2}}\,{\rm d}x ,\end{align}

where now $w=x\eta \lambda$ and (B1) becomes

(B5) \begin{equation} \delta ({\rho },t)=\frac {R f \Delta p}{\mu }\frac {8}{\pi }\int _\rho ^1 \frac {\sqrt {\eta ^2-\rho ^2}}{\eta }\int _0^1 \frac {x e^{-(x\eta \lambda )^2}}{\sqrt {1-x^2}}\,{\rm d}x {\rm d} \eta .\end{equation}

We identify the innermost integral, following (A15), as $F(\eta \lambda )$ , which has the known asymptotic expansion for large values of its argument (A22). We substitute this expansion in (B5) to find the asymptotic expansion of the outer solution for slip, $\delta \sim \delta _{\textit{out}}$ , where

(B6) \begin{equation} \delta _{\textit{out}}(\rho ,t)=\frac {R f \Delta p}{\mu }\frac {8}{\pi } \sum _{n=0}^\infty \frac {1}{\lambda ^{2n+2}} \frac {1}{\sqrt {\pi }}\frac {\varGamma (n+3/2)}{2n+1} \int _\rho ^1 \frac {\sqrt {\eta ^2-\rho ^2}}{\eta ^{2n+3}}\,{\rm d}\eta. \end{equation}

To evaluate the integral, we use the trigonometric variable substitution $\cos \theta =\rho /\eta$ on

(B7) \begin{align} \int _\rho ^1 \frac {\sqrt {\eta ^2-\rho ^2}}{\eta ^{2n+3}}\,{\rm d}\eta =\frac {1}{\rho ^{2n+1}}\int _0^{\textrm{acos}(\rho )}\cos ^{2n}\theta \sin ^2\theta \,{\rm d}\theta \end{align}

and to evaluate this last integral in turn, we will consider the general integral

(B8) \begin{equation} I(m,p,\rho )=\int _0^{\textrm{acos}(\rho )}\cos ^{m}\theta \sin ^p\theta \, {\rm d}\theta \end{equation}

whose solution is as follows. We begin by rearranging the integrand

(B9) \begin{align} I(m,p,\rho )&=\int _0^{\textrm{acos}(\rho )}\left (\cos ^{m-1}\theta \right )\left ( \sin ^p\theta \cos \theta \right ){\rm d}\theta \nonumber \\ &=\frac {\rho ^{m-1}(1-\rho ^2)^{\frac {p+1}{2}}}{p+1}+\frac {m-1}{p-1}\int _0^{\textrm{acos}(\rho )}\cos ^{m-2}\theta \sin ^{p+2}{\rm d}\theta \nonumber \\ &=\frac {\rho ^{m-1}(1-\rho ^2)^{\frac {p+1}{2}}}{p+1}-\frac {m-1}{p+1}\int _0^{\textrm{acos}(\rho )}\cos ^m\theta \sin ^p\theta {\rm d}\theta \nonumber \\ &\quad +\frac {m-1}{p+1}\int _0^{\textrm{acos}(\rho )}\cos ^{m-2}\theta \sin ^p\theta {\rm d}\theta \end{align}

and, in (B9), we pass from the first to second line using integration by parts, and from the second to third via $\sin ^2\theta =1-\cos ^2\theta$ . We identify the first integral as $I(m,p,\rho )$ and the second as $I(m-2,p,\rho )$ . Solving for $I(m,p,\rho )$ , we find the recursion relation:

(B10) \begin{equation} I(m,p,\rho )=\frac {\rho ^{m-1}(1-\rho ^2)^{\frac {p+1}{2}}}{m+p}+\frac {m-1}{m+p}I(m-2,p,\rho ) .\end{equation}

Returning to the specific case $m=2n$ and $p=2$ , we find the recursion relation for which we need only calculate the case $n=0$ to find all others $n\gt 0$ :

(B11) \begin{align} I(0,2,\rho )&=\frac {1}{2}\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right ),\nonumber \\ I(2,2,\rho )&=\frac {1}{8}\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )+\frac {\rho \big(1-\rho ^2 \big)^{3/2}}{4}\nonumber \\ &=\frac {1}{8}\left [\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )+\frac {(1-\rho ^2)^{3/2}}{\rho } 2\rho ^2\right ]\!,\nonumber \\ I(4,2,\rho )&=\frac {1}{16}\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )+\frac {\rho (1-\rho ^2)^{3/2}}{8}+\frac {\rho ^3(1-\rho ^2)^{3/2}}{6}\nonumber \\ &=\frac {1}{16}\left [\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )+\frac {(1-\rho ^2)^{3/2}}{\rho }\left (2\rho ^2+\frac {8}{3}\rho ^4\right )\right ] .\end{align}

From this recursion, we may deduce an explicit expression for $I(m,2,\rho )$ , recognising from (B11) that such an expression will have the form

(B12) \begin{align} I(m,2,\rho )= a_m\left [\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )+\frac { \big(1-\rho ^2 \big)^{3/2}}{\rho }\sum _{k=2,4,\ldots }^mb_k\rho ^k\right ] ,\end{align}

where $a_m$ and $b_k$ are to be determined. From the factor $(m-1)/(m+2)$ in the recursion, we recognise that

(B13) \begin{align} a_m=\frac {(m-1)}{(m+2)}\boldsymbol{\cdot }\frac {(m-3)}{m}\boldsymbol{\cdot }\frac {(m-5)}{(m-2)}\boldsymbol{\cdot }\cdots \boldsymbol{\cdot }\frac {7}{10}\boldsymbol{\cdot }\frac {5}{8}\boldsymbol{\cdot }\frac {3}{6}\boldsymbol{\cdot }\frac {1}{4} \boldsymbol{\cdot }\frac {1}{2}. \end{align}

Substituting $m=2n$ ,

(B14) \begin{align} a_{2n}&=\frac {(2n-1)(2n-3)(2n-5)\ldots }{(2n+2)(2n)(2n-2)\ldots } \\[-12pt] \nonumber \end{align}
(B15) \begin{align} &=\frac {(2n-1)(2n-3)(2n-5)\ldots }{2^{n+1}(n+1)(n)(n-1)\ldots } \\[-12pt] \nonumber \end{align}
(B16) \begin{align} &=\frac {(2n-1)!/[2^{n-1}(n-1)!]}{2^{n+1}(n+1)!} \\[-12pt] \nonumber \end{align}
(B17) \begin{align} &=\frac {(2n-1)!}{2^{2n}(n+1)!(n-1)!}. \end{align}

Passing to the gamma function $n!=\varGamma (n+1)$ ,

(B18) \begin{align} a_{2n} = \frac {\varGamma (2n)}{2^{2n}\varGamma (n)\varGamma (n+2)}. \end{align}

Using the Legendre duplication identity, $\varGamma (2z)/[2^{2z}\varGamma (z)]=\varGamma (z+1/2)/(2\sqrt {\pi })$ ,

(B19) \begin{align} a_{2n}=\frac {1}{2\sqrt {\pi }}\frac {\varGamma (n+1/2)}{\varGamma (n+2)}, \\[-12pt] \nonumber \end{align}
(B20) \begin{align} a_m=\frac {1}{2\sqrt {\pi }}\frac {\varGamma \left (\dfrac {m}{2}+\dfrac {1}{2}\right )}{\varGamma \left (\dfrac {m}{2}+2\right )}. \end{align}

Regarding the coefficients of the series $b_k$ , we can immediately deduce that $b_0=0$ , and may find from the recursion (B10) and the factoring of $a_m$ that, for $k\geqslant 2$ ,

(B21) \begin{align} b_k=\frac {1}{a_k} \frac {1}{k+2}. \end{align}

We may now write the asymptotic expansion for the outer solution under critically stressed conditions:

(B22) \begin{align} \frac {\delta _{\textit{out}}(\rho ,t)}{R f \Delta p/\mu }&=\frac {8}{\pi }\sum _{n=0}^\infty \frac {1}{\lambda ^{2n+2}}\frac {1}{\sqrt {\pi }}\frac {\varGamma (n+3/2)}{2n+1}\frac {1}{\rho ^{2n+1}}a_{2n}\left [\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )\right .\nonumber\\&\quad\left .+\,\frac {(1-\rho ^2)^{3/2}}{\rho }\sum _{k=1}^{n} b_{2k}\rho ^{2k}\right ] ,\end{align}
(B23) \begin{align} \frac {\delta _{\textit{out}}(\rho ,t)}{R f \Delta p/\mu }&=\frac {4}{\pi ^{2}}\sum _{n=0}^\infty \frac {1}{\lambda ^{2n+2}}\frac {1}{\rho ^{2n+1}}\frac {\varGamma (n+3/2)}{2n+1}\frac {\varGamma (n+1/2)}{\varGamma (n+2)} \nonumber\\[3pt] &\quad\boldsymbol{\cdot }\left [\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )+2\sqrt {\pi }\frac {(1-\rho ^2)^{3/2}}{\rho }\sum _{k=1}^{n}\frac {\rho ^{2k}}{2k+2} \frac {\varGamma (k+2)}{\varGamma (k+1/2)}\right ]\!. \end{align}

Appendix C. Alternative derivation of critically stressed outer solution using multipole expansion

We begin by using the inverted solution for slip in terms of the change in shear stress within the crack, repeated as

(C1) \begin{equation} \delta ({\rho },t)=-\frac {4}{\pi }\frac {R}{\mu }\int _\rho ^1 \sqrt {\eta ^2-\rho ^2}\frac {d}{{\rm d}\eta }\int _0^1 \frac {x \Delta \tau (x\eta R,t)}{\sqrt {1-x^2}}\,{\rm d}x{\rm d}\eta ,\end{equation}

where we recall that the change in shear stress from the background level $\Delta \tau (r,t)=\tau _b-\tau (r,t)$ and that we require that the shear stress equals the shear strength $\tau _s=f[\sigma -p(r,t)]$ , with the fluid pressure distribution being

(C2) \begin{equation} p(r,t)=p_o+\Delta pE_1[r^2/(\alpha t)]. \end{equation}

In the critically stressed regime, in which the fluid pressure distribution appears concentrated about the crack origin if viewed over a length scale comparable to the crack size, we approximate the fluid pressure distribution using a multipole expansion

(C3) \begin{equation} p(r,t)=p_o+\Delta p \left [\frac {p_0(t)}{2\pi r}\delta _D(r)-\frac {p_1(t)}{2\pi r}\delta _D^{\prime}(r)+\frac {p_2(t)}{2\pi r}\delta _D^{\prime\prime}(r)-\cdots +\frac {p_n(t)}{2\pi r}\delta _D^{(n)}(r)\right ]\!, \end{equation}

where $\delta _D(r)$ is the Dirac delta, $\delta ^{(n)}_D(r)$ denotes its $n$ th derivative and $p_n(t)$ is the $n$ th moment of the distribution of the scaled, axisymmetric pressure change

(C4) \begin{equation} p_n(t)=\int _0^\infty \frac {r^n}{n!} \frac {p(r,t)-p_o}{\Delta p} 2\pi r {\rm d}r ,\end{equation}

which, for the change due to injection (C2), evaluates to

(C5) \begin{equation} p_n(t)= (\alpha t)^{1+n/2}\,\frac {\varGamma (1+n/2)}{(1+n/2)}\frac {\pi }{n!}. \end{equation}

While $p_n$ will generally be non-zero, we will find that only the even terms in the multipole expansion of $p(r,t)$ will have non-zero contributions to the slip profile. We also recall that the Dirac delta has the property

(C6) \begin{equation} \int _0^\infty \delta ^{(n)}(r a)f(r)\,{\rm d}r=\frac {(-1)^n}{a^{n+1}}f^{(n)}(0). \end{equation}

Expanding the expression (C1),

(C7) \begin{equation} \delta (\rho ,t)= -\frac {4}{\pi }\frac {R}{\mu } \int _\rho ^1 \sqrt {\eta ^2-\rho ^2}\frac {d}{{\rm d}\eta }\int _0^1\frac {x \left [\tau _b-f(\sigma -p(x\eta R,t))\right ] }{\sqrt {1-x^2}}\,{\rm d}x\,{\rm d}\eta, \end{equation}

and rearranging,

(C8) \begin{align} \delta (\rho ,t)&= -\frac {4}{\pi }\frac {R f\Delta p}{\mu } \int _\rho ^1\sqrt {\eta ^2-\rho ^2}\frac {\rm d}{{\rm d}\eta }\int _0^1\frac {x}{\sqrt {1-x^2}} \left [\frac {\tau _b-f\sigma '}{f\Delta p}\right .\nonumber\\&\quad\left . +\, \frac {p(x\eta R,t)-p_o}{\Delta p}\right ] {\rm d}x\,{\rm d}\eta .\end{align}

Considering that a derivative with respect to $\eta$ operates on the innermost integral, we may drop the first term in brackets (identically $T$ ) as it is independent of $\eta$ . Additionally, we can introduce the multipole expansion for the second term of (C3) in brackets to calculate the outer solution asymptotically approximating slip, $\delta \sim \delta _{\textit{out}}$ ,

(C9) \begin{align}& \delta _{\textit{out}}(\rho ,t)= -\frac {4}{\pi }\frac {R f\Delta p}{\mu } \nonumber\\&\quad\boldsymbol{\cdot}\left (\frac {1}{2\pi }\int _\rho ^1\sqrt {\eta ^2-\rho ^2}\frac {\rm d}{{\rm d}\eta }\left [\frac {1}{\eta R}\int _0^1\frac {1}{\sqrt {1-x^2}} \sum _{n=0}^\infty (-1)^n p_n(t)\delta _D^{(n)}(x\eta R) {\rm d}x\right ]{\rm d}\eta \right ). \end{align}

We may use the property of the Dirac delta (C6) to evaluate the innermost integral

(C10) \begin{align}& \delta _{\textit{out}}(\rho ,t)= \frac {4}{\pi }\frac {R f\Delta p}{\mu } \nonumber\\&\quad\boldsymbol{\cdot}\left (\frac {1}{2\pi }\int _\rho ^1\sqrt {\eta ^2-\rho ^2}\frac {\rm d}{{\rm d}\eta }\left [\frac {1}{\eta R}\sum _{n=0}^\infty \frac {p_n(t)}{(\eta R)^{n+1}} \frac {\rm d}{{\rm d}x^n}\left .\left (\frac {1}{\sqrt {1-x^2}}\right )\right |_{x=0}\right ]{\rm d}\eta \right ). \end{align}

We can immediately recognise that odd $n$ th derivatives of the even function $1/\sqrt {1-x^2}$ evaluated at the origin will be zero, due to symmetry. Thus, all odd terms in the sum have no contribution. Rather than continuing the sum as $n=0,2,4\ldots$ , we pass the summation index $n\rightarrow 2n$ with the new index maintaining the progression $n=0,1,2,\ldots$ . From the Maclaurin series (A16), we find that, for this new index $n$ ,

(C11) \begin{align} \frac {\rm d}{{\rm d}x^n}\left .\left (\frac {1}{\sqrt {1-x^2}}\right )\right |_{x=0}=\left (\frac {(2n)!}{2^{n} \,n !}\right )^2, \end{align}

such that

(C12) \begin{equation} \delta _{\textit{out}}(\rho ,t)= \frac {4}{\pi }\frac {R f\Delta p}{\mu } \left (\frac {1}{2\pi }\int _\rho ^1\sqrt {\eta ^2-\rho ^2}\frac {\rm d}{{\rm d}\eta }\left [\frac {1}{\eta R}\sum _{n=0}^\infty \frac {p_{2n}(t)}{(\eta R)^{2n+1}} \left (\frac {(2n)!}{2^{n} \,n!}\right )^2\right ]{\rm d}\eta \right ) .\end{equation}

Evaluating the derivative, we arrive to

(C13) \begin{equation} \delta _{\textit{out}}(\rho ,t)= \frac {4}{\pi }\frac {R f\Delta p}{\mu } \left (\frac {1}{2\pi }\sum _{n=0}^\infty \frac {p_{2n}(t)}{R^{2n+2}}\left (\frac {(2n)!}{2^{n}\,n!}\right )^2 (2n+2)\int _\rho ^1\frac {\sqrt {\eta ^2-\rho ^2}}{\eta ^{2n+3}}\,{\rm d}\eta \right ). \end{equation}

Using (C5) to evaluate $p_{2n}$ and recognising that $(\alpha t)^{1+n}/R^{2n +2}=1/\lambda ^{2n+2}$ , this further simplifies to

(C14) \begin{equation} \delta _{\textit{out}}(\rho ,t)= \frac {4}{\pi }\frac {R f\Delta p}{\mu } \sum _{n=0}^\infty \frac {1}{\lambda ^{2n+2}}\frac {(2n)!}{2^{2n}\,n!} \int _\rho ^1\frac {\sqrt {\eta ^2-\rho ^2}}{\eta ^{2n+3}}\,{\rm d}\eta \end{equation}

with which we may compare the solution for $\delta _{\textit{out}}$ found by direct asymptotic expansion, (B6). The equivalence of these two expressions for $\delta _{\textit{out}}$ follows from the equivalence of the numerical coefficients

(C15) \begin{align} \frac {8}{\pi } \frac {1}{\sqrt {\pi }}\frac {\varGamma (n+3/2)}{2n+1}=\frac {4}{\pi }\frac {(2n)!}{2^{2n} n!} ,\end{align}

which directly follows from the Legendre-duplication formula, (A21). The final expression for the outer solution can then be found by repeating the steps subsequent to (B6).

Appendix D. Asymptotic solution details for marginally pressurised conditions

We pursue the regular perturbation expansion of the solution in the marginally pressurised limit, $\lambda \ll 1$ . We return to the integral expression for slip due to the constant-injection-rate source

(D1) \begin{align} \delta ({\rho },t)=\frac {8}{\pi }\frac {R f\Delta p}{\mu }\int _\rho ^1 \frac {\sqrt {\eta ^2-\rho ^2}}{\eta }\int _0^1 \frac {x e^{-(x\eta \lambda )^2}}{\sqrt {1-x^2}}\,{\rm d}x{\rm d}\eta, \end{align}

where we again identify the innermost integral as $F(\eta \lambda )$ and look to expand it now for $\eta \lambda \ll 1$ . In (D1), we can substitute the Taylor expansion

(D2) \begin{align} e^{-(x\eta \lambda )^2}=\sum _{n=0}^\infty \frac {(-1)^{n}}{n!} (x\eta \lambda )^{2n} \end{align}

and subsequently evaluate integrals of the form

(D3) \begin{align} I_{2n+1}=\int _0^1 \frac {x^{2n+1}}{\sqrt {1-x^2}}\,{\rm d}x \end{align}

for $n=0,1,2\ldots$ Letting $m=2n+1$ and $x=\sin \theta$ , we transform the integral to

(D4) \begin{align} I_m&=\int _0^1 \frac {x^m}{\sqrt {1-x^2}}\,{\rm d}x=\int _0^{\pi /2}\sin ^m\theta \,{\rm d}\theta \nonumber \\ &=\int _0^{\pi /2}\!\sin ^{m-1}\theta \sin \theta \,{\rm d}\theta =\left (-\sin ^{m-1}\theta \cos \theta \right | _0^{\pi /2}+(m-1)\int _0^{\pi /2}\!\sin ^{m-2}\theta \cos ^2\theta {\rm d}\theta \nonumber \\ &=(m-1)\left (\int _0^{\pi /2}\sin ^{m-2}\theta \,{\rm d}\theta -\int _0^{\pi /2}\sin ^{m}\theta \,{\rm d} \theta \right )\nonumber \\ &=(m-1)(I_{m-2}-I_m) \end{align}

and arrive to the recursion relation

(D5) \begin{align} I_m=\frac {m-1}{m}I_{m-2}. \end{align}

Given that $I_1=1$ ,

(D6) \begin{equation} I_{2n+1}=\frac {(2n)\boldsymbol{\cdot }(2n-2)\boldsymbol{\cdot }\cdots \boldsymbol{\cdot }6\boldsymbol{\cdot }4\boldsymbol{\cdot }2}{(2n+1)\boldsymbol{\cdot }(2n-1)\boldsymbol{\cdot }\cdots \boldsymbol{\cdot }7\boldsymbol{\cdot }5\boldsymbol{\cdot }3}= \frac {2^n n!}{(2n+1)!/(2^n n!)}, \end{equation}

which we may rewrite using the Gamma function

(D7) \begin{align} I_{2n+1}=4^n \frac {\left [\varGamma (n+1)\right ]^2}{\varGamma (2n+2)} \end{align}

and using the Legendre duplication formula, we can further simplify

(D8) \begin{align} I_{2n+1}=\frac {\sqrt {\pi }}{2} \frac {\varGamma (n+1)}{\varGamma (n+3/2)} \end{align}

and arrive to an expression for the expansion of $F$ for small values of its argument

(D9) \begin{align}&\qquad\qquad\qquad F(\eta \lambda )= \frac {\sqrt {\pi }}{2}\sum _{n=0}^\infty (\eta \lambda )^{2n} \frac {(-1)^n}{\varGamma (n+3/2)}, \end{align}
(D10) \begin{align} \frac {\delta (\rho ,t)}{R f \Delta p/\mu }&=\frac {8}{\pi }\int _\rho ^1 \!\frac {\sqrt {\eta ^2-\rho ^2}}{\eta } F(\eta \lambda ){\rm d}\eta =\frac {4}{\sqrt {\pi }}\sum _{n=0}^\infty \frac {(-1)^n}{\varGamma (n+3/2)}\lambda ^{2n} \!\int _\rho ^1 \!\sqrt {\eta ^2-\rho ^2}\eta ^{2n-1}{\rm d}\eta. \end{align}

Recycling the symbol $I$ , we denote the integrals within the sum as

(D11) \begin{equation} I_{2n-1}(\rho )=\int _\rho ^1 \sqrt {\eta ^2-\rho ^2}\eta ^{2n-1}\,{\rm d}\eta \end{equation}

for $n=0,1,2\ldots$

For $n=0$ , we may evaluate $I_{-1}$ using the substitution $\rho /\eta =\cos \theta$ and integration by parts

(D12) \begin{align} I_{-1}(\rho )=\int _\rho ^1 \frac {\sqrt {\eta ^2-\rho ^2}}{\eta }{\rm d}\eta &=\rho \int _0^{\textit{acos}(\rho )} \frac {\sin ^2\theta } {\cos ^2\theta }\,{\rm d}\theta \notag \nonumber \\ &=\rho \left [\left (\frac {\sin \theta }{\cos \theta }\right | _0^{\textit{acos}(\rho )}-\int _0^{\textit{acos}(\rho )}\,{\rm d}\theta \right ]\notag \nonumber \\ &=\sqrt {1-\rho ^2}-\rho \,\textrm{acos}(\rho ) .\end{align}

The remaining integrals in the sum ( $n\gt 0$ ) can be found by recursion. Letting $m=2n-1$ and performing integration by parts,

(D13) \begin{align} I_{m}(\rho )&=\int _\rho ^1 \sqrt {\eta ^2-\rho ^2}\eta ^{m}\,{\rm d}\eta \nonumber \\&=\frac {1}{3}\left . \big(\eta ^2-\rho ^2 \big)^{3/2}\eta ^{m-1}\right |_\rho ^1-\frac {m-1}{3}\int _\rho ^1 \big(\eta ^2-\rho ^2 \big)^{3/2}\eta ^{m-2}\,{\rm d}\eta \nonumber \\&=\frac {1}{3}\big(1-\rho ^2 \big)^{3/2}-\frac {m-1}{3}\int _\rho ^1 \big(\eta ^2-\rho ^2 \big)^{1/2}\eta ^m \big(\eta ^2-\rho ^2 \big)\eta ^{-2}\,{\rm d}\eta \nonumber \\&=\frac {1}{3}(1-\rho ^2)^{3/2}-\frac {m-1}{3}\big (I_m(\rho )-\rho ^2 I_{m-2}(\rho )\big ) ,\end{align}

which leads to the recursion relation

(D14) \begin{align} I_{m}(\rho )=\frac {\big(1-\rho ^2 \big)^{3/2}}{m+2}+\frac {m-1}{m+2}\rho ^2 I_{m-2}(\rho ) \end{align}

with $I_1(\rho )=(1-\rho ^2)^{3/2}/3$ . We may solve this recurrence relation to yield an explicit expression for $I_m(\rho )$ . Given the recursion relation, we can expect a solution of the form

(D15) \begin{align} I_m(\rho )=\frac {\big(1-\rho ^2 \big)^{3/2}}{m+2}\sum _{k=1,3,5\ldots }^m a_{km}\rho ^{k-1}, \end{align}

where

(D16) \begin{align} &a_{1m}=1,\nonumber \\ &a_{3m}=\frac {(m-1)}{m},\nonumber \\ &a_{5m}=\frac {(m-1)}{m}\boldsymbol{\cdot }\frac {(m-3)}{(m-2)},\nonumber \\& a_{7m}=\frac {(m-1)}{m}\boldsymbol{\cdot }\frac {(m-3)}{(m-2)}\boldsymbol{\cdot }\frac {(m-5)}{(m-4)},\nonumber \\ &\ldots \nonumber \\ &a_{mm}=\frac {(m-1)}{m}\boldsymbol{\cdot }\frac {(m-3)}{(m-2)}\boldsymbol{\cdot }\frac {(m-5)}{(m-4)}\boldsymbol{\cdot }\cdots \boldsymbol{\cdot }\frac {2}{3}. \end{align}

For convenience, we pass from odd-numbered $m$ and $k$ to a regular series of integers via $m=2n+1$ and $k=2j+1$ with $n$ and $j$ taking on the values $0,1,2,3,\ldots$ such that we may write the integral $I_{m}$ as

(D17) \begin{align} I_{2n+1}(\rho )=\frac {(1-\rho ^2)^{3/2}}{2n+3}\sum _{j=0}^n c_{jn}\rho ^{2j} \end{align}

with

(D18) \begin{align} c_{0n}&=1,\nonumber \\ c_{1n}&=\frac {(2n)}{(2n+1)},\nonumber \\ c_{2n}&=\frac {(2n)}{(2n+1)}\boldsymbol{\cdot }\frac {(2n-2)}{(2n-1)},\nonumber \\ c_{3n}&=\frac {(2n)}{(2n+1)}\boldsymbol{\cdot }\frac {(2n-2)}{(2n-1)}\boldsymbol{\cdot }\frac {(2n-4)}{(2n-3)},\nonumber \\ \ldots \nonumber \\ c_{nn}&=\frac {(2n)}{(2n+1)}\boldsymbol{\cdot }\frac {(2n-2)}{(2n-1)}\boldsymbol{\cdot }\frac {(2n-4)}{(2n-3)}\boldsymbol{\cdot }\cdots \boldsymbol{\cdot }\frac {2}{3}. \end{align}

Given experience with similar products in preceding appendices, we may express this as

(D19) \begin{equation} c_{jn}=\frac {2^n n!}{2^{n-j}(n-j)!} \boldsymbol{\cdot }\frac {(2(n-j)+1)!/(2^{n-j}(n-j)!)}{(2n+1)!/(2^n n!)}, \end{equation}

where the term before $\boldsymbol{\cdot}$ provides the products found in the numerator of the preceding expressions for $c_j$ , whereas the term following $\boldsymbol{\cdot}$ provides that found in the denominator. Proceeding to further simplify this expression for $c_j$ ,

(D20) \begin{align} c_{jn}&=\frac {2^{2n} (n!)^2}{2^{2(n-j)}((n-j)!)^2} \boldsymbol{\cdot }\frac {(2(n-j)+1)!}{(2n+1)!}\nonumber \\ &=\frac {\varGamma (n+1)}{\varGamma (n+3/2)}\frac {\varGamma (n-j+3/2)}{\varGamma (n-j+1)}, \end{align}

where the first line simply rearranges (D19) and the latter makes use of the Gamma function and its the Legendre duplication formula (A21). We thus arrive at the sought-after expression for $I_{2n+1}$ , for $n=0,1,2\ldots\,$ . However, to be consistent with our use of $n$ in our original definition of the integral $I$ , (D11), we shift $n$ and $j$ by one, such that $n$ and $j$ now take on values $1,2,3,\ldots\,$ , and

(D21) \begin{align} I_{2n-1}(\rho )=\frac {(1-\rho ^2)^{3/2}}{2n+3} \frac {\varGamma (n)}{\varGamma (n+1/2)}\sum _{j=1}^n \frac {\varGamma (n-j+3/2)}{\varGamma (n-j+1)}\rho ^{2j-2}. \end{align}

Such that, upon substitution of this expression and that for $I_{-1}(\rho )$ , (D12), into the expression for slip in the marginally pressurised limit, $\lambda \rightarrow 0$ ,

(D22) \begin{align} \frac {\delta (\rho ,t)}{R(t)f \Delta p/\mu }=\frac {8}{\pi }\int _{\rho }^1 \frac {\sqrt {\eta ^2-\rho ^2}}{\eta }{\rm d}\eta +\frac {4}{\sqrt {\pi }}\sum _{n=1}^\infty \lambda ^{2n}\frac {(-1)^n}{\varGamma (n+3/2)}\int _{\rho }^1 \sqrt {\eta ^2-\rho ^2}\eta ^{2n-1}{\rm d}\eta \end{align}

leads to

(D23) \begin{align} &\frac {\delta (\rho ,t)}{R(t)f \Delta p/\mu }=\frac {8}{\pi }\left (\sqrt {1-\rho ^2}-\rho \text{ acos}(\rho )\right )\nonumber \\&\quad+\frac {4}{\sqrt {\pi }}\big (1-\rho ^2\big )^{3/2}\sum _{n=1}^\infty \frac {\lambda ^{2n}}{2n+1}\frac {(-1)^n}{\varGamma (n+3/2)}\frac {\varGamma (n)}{\varGamma (n+1/2)}\sum _{j=1}^n \frac {\varGamma (n-j+3/2)}{\varGamma (n-j+1)}\rho ^{2j-2}. \end{align}

Appendix E. Far-field expansion of inner solution

In the interest of determining the overlap of the inner and outer solutions, here, we derive the asymptotic expansion of the inner solution (3.1)–(3.2) in the far field ( $s\rightarrow \infty$ ). Specifically, we focus on the leading-order solution (3.2), reproduced as

(E1) \begin{align} \delta _0(s,t)=2\sqrt {\pi }e^{-s^2/2} \big [\big(1+s^2 \big)I_0 \big(s^2/2 \big)+s^2 I_1 \big(s^2/2 \big)\big ]-4s, \end{align}

which follows from its integral representation (A11)

(E2) \begin{equation} \delta _0(s,t)=\frac {8}{\pi }\int _s^\infty \sqrt {1-(s/\xi )^2}F(\xi )\,{\rm d}\xi, \end{equation}

where $F(\xi )$ is given by (A7).

We begin by substituting the asymptotic expansion of $F(\xi )$ for large values of its argument, (A22),

(E3) \begin{equation} F(\xi )\sim \sum _{n=0}^\infty \frac {1}{\sqrt {\pi }}\frac {1}{\xi ^{2n+2}}\frac {\varGamma (n+3/2)}{2n+1}, \end{equation}

into (E2), which then requires us to evaluate integrals of the form

(E4) \begin{align} \int _s^\infty \frac {\sqrt {1-(s/\xi )^2}}{\xi ^{2n+2}}\,{\rm d}\xi, \end{align}

which suggests the change of variable $\cos \theta =s/\xi$ , transforming this integral to

(E5) \begin{align} \frac {1}{s^{2n+1}}\int _0^{\pi /2}\sin ^2\theta \cos ^{2n}\theta \,{\rm d}\theta, \end{align}

the integral that is a special case of the one encountered in Appendix B,

(E6) \begin{align} I(m,p,\rho )=\int _0^{\textrm{acos}(\rho )}\cos ^m\theta \sin ^p\theta \, {\rm d}\theta. \end{align}

The particular solution for $p=2$ was found there to be

(E7) \begin{align} I(m,2,\rho )&= \frac {1}{2\sqrt {\pi }}\frac {\varGamma \!\left (\dfrac {m}{2}+\dfrac {1}{2}\right )}{\varGamma \!\left (\dfrac {m}{2}+2\right )}\left [\left (\textrm{acos}(\rho )-\rho \sqrt {1-\rho ^2}\right )\vphantom{\frac {\varGamma \!\left (\dfrac {k}{2}+2\right )}{\varGamma \!\left (\dfrac {k}{2}+\dfrac {1}{2}\right )}}\right .\nonumber\\&\quad\left .+\,\frac {(1-\rho ^2)^{3/2}}{\rho }\sum _{k=2,4,\ldots }^m \frac {2\sqrt {\pi }}{k+2}\frac {\varGamma \!\left (\dfrac {k}{2}+2\right )}{\varGamma \!\left (\dfrac {k}{2}+\dfrac {1}{2}\right )}\rho ^k\right ]\!, \end{align}

which simplifies for $\rho =0$ and $m=2n$ to

(E8) \begin{align} I(2n,2,0)=\frac {\sqrt {\pi }}{4}\frac {\varGamma (n+1/2)}{\varGamma (n+2)} .\end{align}

Thus, following the substitution of (E3) into (E2), we deduce the asymptotic expansion for $\delta _0$ as $s\rightarrow \infty$ ,

(E9) \begin{align} \delta _0(s,t)&\sim \frac {8}{\pi }\sum _{n=0}^\infty \frac {1}{\sqrt {\pi }}\frac {\varGamma (n+3/2)}{2n+1}\int _s^\infty \frac {\sqrt {1-(s/\xi )^2}}{\xi ^{2n+2}}\,{\rm d}\xi \\[3pt] &\sim \frac {2}{\pi }\sum _{n=0}^\infty \frac {\varGamma (n+3/2)}{2n+1}\frac {\varGamma (n+1/2)}{\varGamma (n+2)}\frac {1}{s^{2n+1}}. \end{align}

References

Atkinson, C. & Craster, R.V. 1994 A singular perturbation approach to integral equations occurring in poroelasticity. IMA J. Appl. Maths 52, 221–252.CrossRefGoogle Scholar
Bender, C.M. & Orszag, S.A. 1999 Advanced mathematical methods for scientists and engineers. Springer-Verlag.10.1007/978-1-4757-3069-2CrossRefGoogle Scholar
Berrios-Rivera, N., Ozawa, S. & Dunham, E.M. 2025 Models of fluid-induced seismic slip with permeability enhancement and rate-and-state friction. ESS Open Archive. Available at: https://essopenarchive.org/doi/full/10.22541/essoar.175357075.53109087/v1.Google Scholar
Bhattacharya, P. & Viesca, R.C. 2019 Fluid-induced aseismic fault slip outpaces pore-fluid migration. Science 364, 464–468.CrossRefGoogle ScholarPubMed
Craster, R.V. & Atkinson, C. 1992 Shear cracks in thermoplastic and poroelastic media. J. Mech. Phys. Solids 40 (4), 887–924.CrossRefGoogle Scholar
Craster, R.V. & Atkinson, C. 1994 Crack problems in a poroelastic medium: an asymptotic approach. Proc. R. Soc. Lond. A 346, 387–428.Google Scholar
Desroches, J., Detournay, E., Lenoach, B., Papanastasiou, P., Pearson, J.R.A., Thiercelin, M. & Cheng, A. 1994 The crack tip region in hydraulic fracturing. Proc. R. Soc. Lond. A 447, 39–48.Google Scholar
Dublanchet, P. 2019 Fluid driven shear cracks on a strengthening rate-and-state frictional fault. J. Mech. Phys. Solids 132 , 103672.CrossRefGoogle Scholar
Dunham, E.M. 2024 Fluid-driven aseismic fault slip with permeability enhancement and dilatancy. Phil. Trans. R. Soc. A: Math. Phys. Eng. Sci. 382, 20230255.Google ScholarPubMed
Dyskin, A.V., Germanovich, L.N. & Ustinov, K.B. 2000 Asymptotic analysis of crack interaction with free boundary. Intl J. Solids Struct 37, 857–886.10.1016/S0020-7683(99)00063-3CrossRefGoogle Scholar
Edmunds, T.M. & Willis, J.R. 1976 a Matched asymptotic expansions in nonlinear fracture mechanics – I. Longitudinal shear of an elastic perfectly-plastic specimen. J. Mech. Phys. Solids 24 (4), 205–223.Google Scholar
Edmunds, T.M. & Willis, J.R. 1976 b Matched asymptotic expansions in nonlinear fracture mechanics – II. Longitudinal shear of an elastic work-hardening plastic specimen. J. Mech. Phys. Solids 24 (4), 225–237.10.1016/0022-5096(76)90004-1CrossRefGoogle Scholar
Edmunds, T.M. & Willis, J.R. 1977 Matched asymptotic expansions in nonlinear fracture mechanics – III. In-plane loading of an elastic perfectly-plastic symmetric specimen. J. Mech. Phys. Solids 25 (6), 423–455.CrossRefGoogle Scholar
Erickson, B.A. et al. 2020 The community code verification exercise for simulating sequences of earthquakes and aseismic slip (SEAS). Seismol. Res. Lett. 91 (2A), 874–890.10.1785/0220190248CrossRefGoogle Scholar
Garagash, D.I. 2006 Propagation of a plane-strain hydraulic fracture with a fluid lag: early-time solution. Intl J. Solids Struct. 43, 5811–5835.10.1016/j.ijsolstr.2005.10.009CrossRefGoogle Scholar
Garagash, D.I. & Detournay, E. 2000 The tip region of a fluid-driven fracture in an elastic medium. J. Appl. Mech. 67, 183–192.Google Scholar
Garagash, D.I. & Detournay, E. 2005 Plane-strain propagation of a fluid-driven fracture: small toughness solution. J. Appl. Mech 72 (6), 916–928.Google Scholar
Garagash, D.I. & Germanovich, L.N. 2012 Nucleation and arrest of dynamic slip on a pressurized fault. J. Geophys. Res. 117 , B10310.Google Scholar
Harris, R.A. et al. 2009 The SCEC/USGS dynamic earthquake rupture code verification exercise. Seismol. Res. Lett. 80 (1), 119–126.CrossRefGoogle Scholar
Hinch, E.J. 1991 Perturbation Methods. Cambridge University Press.CrossRefGoogle Scholar
Im, K. & Avouac, J.-P. 2024 Maximum magnitude of induced earthquakes in rate and state friction framework. Seismol. Res. Lett. 96 (3), 1654–1664.10.1785/0220240382CrossRefGoogle Scholar
Jacquey, A.B. & Viesca, R.C. 2023 Nucleation and arrest of fluid-induced aseismic slip. Geophys. Res. Lett. 50 , e2022GL101228.CrossRefGoogle Scholar
Lambert, V.R. et al. 2025 Community-driven code comparisons for simulations of fluid-induced aseismic slip. J. Geophys. Res 130, e2024JB030601.10.1029/2024JB030601CrossRefGoogle Scholar
Lister, J.R. 1990 Buoyancy-driven fluid fracture: similarity solutions for the horizontal and vertical propagation of fluid-filled cracks. J. Fluid Mech. 217, 213–239.10.1017/S0022112090000696CrossRefGoogle Scholar
Lister, J.R., Skinner, D.J. & Large, T.M.J. 2019 Viscous control of shallow elastic fracture: peeling without precursors. J. Fluid. Mech 868, 119–140.10.1017/jfm.2019.185CrossRefGoogle Scholar
Marck, J., Savitski, A.A. & Detournay, E. 2015 Line source in a poroelasic layer bounded by an elastic space. Intl J. Numer. Anal. Meth. Geomech. 39, 1484–1505.10.1002/nag.2405CrossRefGoogle Scholar
Rice, J.R. 1980 Mechanics of earthquake rupture. In Physics of the Earth’s Interior (eds. A.M. Dziewonski & E. Boschi), pp. 555–649. Italian Physical Society/North Holland Publ. Co.Google Scholar
Ross, Z.E., Cochran, E.S., Trugman, D.T. & Smith, J.D. 2020 3D fault architecture controls the dynamism of earthquake swarms. Science 368, 1357–1361.10.1126/science.abb0779CrossRefGoogle ScholarPubMed
Rudnicki, J.W. 1990 Boundary layer analysis of plane strain shear cracks propagating steadily on an impermeable plane in an elastic diffusive solid. J. Mech. Phys. Solids 39 (2), 201–221.10.1016/0022-5096(91)90003-7CrossRefGoogle Scholar
Rudnicki, J.W. & Koutsibelas, D.A. 1991 Steady propagation of plane strain shear cracks on an impermeable plane in an elastic diffusive solid. Intl J. Solids Struct. 27 (2), 205–225.CrossRefGoogle Scholar
Sáez, A. & Lecampion, B. 2023 Post-injection aseismic slip as a mechanism for the delayed triggering of seismicity. Proc. R. Soc. A: Math. Phys. Eng. Sci. 479, 20220810.10.1098/rspa.2022.0810CrossRefGoogle Scholar
Sáez, A. & Lecampion, B. 2024 Fluid-Driven slow slip and earthquake nucleation on a slip-weakening circular fault. J. Mech. Phys. Solids 183, 105506. 10.1016/j.jmps.2023.105506CrossRefGoogle Scholar
Sáez, A., Lecampion, B., Bhattacharya, P. & Viesca, R.C. 2022 Three-dimensional fluid-driven stable frictional ruptures. J. Mech. Phys. Solids 160, 104754.10.1016/j.jmps.2021.104754CrossRefGoogle Scholar
Sáez, A., Passelègue, F. & Lecampion, B. 2025 Maximum size and magnitude of injection-induced slow slip events. Sci. Adv. 11, eadq0662.10.1126/sciadv.adq0662CrossRefGoogle ScholarPubMed
Salamon, N.J. & Dundurs, J. 1971 Elastic fields of a dislocation loop in a two-phase material. J. Elasticity 1 (2), 153–164.Google Scholar
Salamon, N.J. & Dundurs, J. 1977 A circular glide dislocation loop in a two-phase material. J. Phys. C Solid State 10, 497–507.10.1088/0022-3719/10/4/007CrossRefGoogle Scholar
Savitski, A.A. & Detournay, E. 2002 Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions. Intl J. Solids Struct. 39, 6311–6337.10.1016/S0020-7683(02)00492-4CrossRefGoogle Scholar
Simons, D.A. 1977 Boundary-layer analysis of propagating mode II cracks in porous elastic media. J. Mech. Phys. Solids 25, 99–115.CrossRefGoogle Scholar
Simons, D.A. & Rice, J.R. 1976 The stabilization of spreading shear faults by coupled deformation-diffusion effects in fluid-infiltrated porous materials. J. Geophys. Res. 81 (29), 5322–5334.Google Scholar
Sneddon, I.N. 1951 Fourier Transforms. McGraw-Hill Google Scholar
Spence, D.A. & Sharp, P. 1985 Self-similar solutions for elastohydrodynamic cavity flow. Proc. R. Soc. Lond. A 400, 289–313.Google Scholar
Thouless, M.D., Evans, A.G., Ashby, M.F. & Hutchinson, J.W. 1987 The edge cracking and spalling of brittle plates. Acta Metall. 35 (6), 1333–1341.10.1016/0001-6160(87)90015-0CrossRefGoogle Scholar
Viesca, R.C. 2021 Self-similar fault slip in response to fluid injection. J. Fluid Mech. 928, A29.10.1017/jfm.2021.825CrossRefGoogle Scholar
Viesca, R.C. & Rice, J.R. 2012 Nucleation of slip-weakening rupture instability in landslides by localized increase of pore pressure. J. Geophys. Res. 117, B03104.Google Scholar
Zlatin, A.N. & Khrapkov, A.A. 1986 A semi-infinite crack parallel to the boundary of the elastic half-plane. Sov. Dokl. Phys. 31 (12), 1009–1010.Google Scholar
Figure 0

Figure 1.Plot of the variation of the stress-injection parameter $T$ (solid black), which characterises the pre-injection state of fault stress and strength as well as the strength of the injection, with the prefactor $\lambda$ for the self-similar expansion of the crack radius, $R(t)=\lambda \sqrt {\alpha t}$. The marginally pressurised (2.21, blue-dashed) and critically stressed (2.19, red-dashed) are the perturbation expansions in the small and large $\lambda$ limit, respectively. The marginally pressurised expansion is carried to $O(\lambda ^8)$ and the critically stressed expansion is carried to $O(\lambda ^{-4})$, which represents an optimal truncation for $\lambda$ near 1 (discussed further in § 6.1).

Figure 1

Figure 2.Plots of inner (red dashed) and outer (blue dashed) solutions in comparison with full solution (solid black), for a case modestly within the critically stressed regime, $\lambda =2$. The inner and outer solutions are carried out to (a) $O(\lambda ^{-1})$, (b) $O(\lambda ^{-3})$, (c) $O(\lambda ^{-5})$, (d) $O(\lambda ^{-7})$ and (e) $O(\lambda ^{-9})$. Panel (e) shows the optimally truncated inner and outer solutions of the expansion (discussion further in § 6.1).

Figure 2

Figure 3.Same as in figure 2, except with the ordinate axes now on a logarithmic scale to highlight separation of inner and outer solutions, as well as their diminishing overlap as order is increased to the order of optimal truncation (panel e).

Figure 3

Figure 4.(a) Comparison of composite solutions for slip at various orders $\delta _{cs}$ (dashed lines) against the solution for slip $\delta$ (black line) at a value of $\lambda =2$ only modestly in the critically stressed regime. Composite solutions are carried out to the following orders in the small parameter $\lambda ^{-1}$: leading first (red), third (blue), fifth (green), seventh (cyan) and ninth (orange). By fifth order, the composite solution is indistinguishable from the true solution. (b) Absolute error of composite solutions at various orders are shown as a solid line, whose colours indicate order of composite solution and correspond to those of dashed lines in panel (a). Dashed and dotted lines indicate value of the next-order term in the composite solution: e.g. the red-dashed line shows that the value of the third-order ($\lambda ^{-3}$) term of the composite solution is equal to more than 50 % of the error of the first-order solution (red-solid line). Note the dramatic reduction in error for the composite solution at ninth order (solid orange line), reflecting optimal truncation of the composite solution (discussed further in § 6.1). (c) Deviation of the composite solution $\delta _{cs}$ at various orders (dashed lines) from the solution for slip $\delta$ near the crack centre, $r=0$. Excellent agreement of optimally truncated composite solution at ninth order is again apparent. (d) Deviation of the composite solution near the crack tip, $r=R(t)$.

Figure 4

Figure 5.(a) Comparison of the marginally pressurised series solution truncated at various orders (dashed lines) with the full solution (black solid line), both for $\lambda =1/2$. The series are carried up to: $O(\lambda ^0)$ (red), $O(\lambda ^2)$ (blue), $O(\lambda ^4)$ (green), $O(\lambda ^6)$ (cyan) and $O(\lambda ^8)$ (orange). Beyond the leading order, the marginally pressurised solutions are indistinguishable at this scale. (b) Plot of the relative error of the truncated series solutions, with the same colour scale as panel (a). The regular perturbation series exhibits uniform convergence in terms of relative error.

Figure 5

Figure 6.Comparison of (a) aseismic moment and (b) peak slip as a function of $\lambda$ (solid black) compared with truncated marginally pressurised (blue dashed) and critically stressed (red dashed) expansions. (a) The marginally pressurised expansion is $O(\lambda ^{8})$ and the critically stressed expansion is $O(\lambda ^{-6})$. The former reflects the first four terms of (5.2); the latter reflects the first three terms of (5.3) and a $\lambda ^{-1}$ scaling factor due to a change in normalisation. (b) The marginally pressurised expansion $\delta _{\textit{mp}}$ is $O(\lambda ^9)$ and the critically stressed expansion $\delta _{cs}$ is $O(\lambda ^{-5})$. The former reflects the first five terms of (5.5) and a $\lambda$ scaling factor due to a change in normalisation. The latter reflects the first three terms of (5.4). In both panels, truncation of the critically stressed expansions at three terms is optimal for values of $\lambda$ close to 1. Optimal truncation is discussed in § 6.1.

Figure 6

Figure 7.(a) Plot of absolute error at the crack centre of the asymptotic series solution (5.4) truncated at $j=N$, for $\lambda =2$. The optimal truncation is expected to be $N=\lambda ^2$, and we observe a dramatic drop in error for $N=4$. (b) Absolute error of the composite solution at the crack centre optimally truncated at $N=\lambda ^2$, plotted for fixed values of $\lambda$ (black dots). The expected dependence of the error on $\lambda$ (6.6) is shown as a solid line. (Bottom) Spatial distribution of absolute error showing uniform convergence for five composite solutions (3.9) optimally truncated at $N=\lambda ^2$: $\lambda =1$ (red), $\lambda =2$ (blue), $\lambda =4$ (green), $\lambda =6$ (cyan) and $\lambda =8$ (orange).

Figure 7

Figure 8.(a) Absolute error of the slip at the crack centre for the marginally pressurised series solution at $O(\lambda ^{2n})$ for the particular case $\lambda =1/2$ (black dots). The expected error (6.10) is shown as a solid line. (b) Distribution of the absolute error along the expanding crack for the case $\lambda =1/2$ truncated at $O(\lambda ^0)$ (red), $O(\lambda ^2)$ (blue), $O(\lambda ^4)$ (green), $O(\lambda ^6)$ (cyan), $O(\lambda ^8)$ (orange).

Figure 8

Figure 9.Plot of the variation of the prefactor $\lambda$ for the self-similar expansion of the crack radius, $R(t)=\lambda \sqrt {\alpha t}$ (solid black), as it depends on the stress-injection parameter $T$, which characterises the pre-injection state of fault stress and strength as well as the strength of the injection. This plot is the inverse of figure 1 and showcases the asymptotic expansions of $\lambda$ in terms of $T$: in comparison are the critically stressed (7.3, red dashed) and the marginally pressurised (7.6, blue dashed) perturbation expansions in the small and large $T$ limit, respectively. The marginally pressurised expansion shown here is carried to $O(e^{(-7/2)T})$ and the critically stressed expansion is carried to $O(T^{(1/2)})$, which represents an optimal truncation for $T$ near 1.