Induced aseismic slip and the onset of seismicity in displaced faults

Abstract

We address aseismic fault slip and the onset of seismicity resulting from depletion-induced or injection-induced stresses in reservoirs with pre-existing vertical or inclined faults. Building on classic results, we discuss semi-analytical modelling techniques for fault slip including dislocation theory, Cauchy-type singular integral equations and the use of Chebyshev polynomials for their solution and an eigenvalue-based stability analysis. We consider slip patch development during depletion for faults with zero, constant static and slip-weakening friction, and our results confirm earlier findings based on numerical simulation, in particular the aseismic growth of two slip patches that may subsequently merge and/or become unstable resulting in nucleation of seismic slip. New findings include improved approximate expressions for the induced seismic moment per unit strike length and a description of the effect of coupling between the slip patches which affects both forward simulation and eigenvalue computation for high values of the ratio of fault throw to reservoir height. Our implementation based on analytical inversion and semi-analytical integration with Chebyshev polynomials is more efficient and more robust than our numerical integration approach. It is not yet well suited for Monte Carlo simulation, which typically requires sub-second simulation times, but with some further development that option seems to be within reach. Moreover, our results offer a possibility for embedded fault modelling in large-scale numerical simulation tools.

Introduction

The motivation for our paper is to contribute to the understanding of earthquakes that result from the production of hydrocarbons, geothermal operations or the storage of fluids in subsurface reservoirs (Segall, 1989; Suckale, 2009; Elsworth et al., 2016; Muntendam-Bos et al., 2021). Therefore, we investigate the use of (semi-) analytical techniques to describe induced fault slip and their potential for computational purposes. Our paper forms an extension to an earlier one that provided closed-form analytical expressions for the injection-induced or depletion-induced elastic stresses in and around displaced faults, that is, faults with a non-zero offset (Jansen et al., 2019). That earlier paper build methodologically on analytical work of, amongst others, Geertsma (1973), Segall (1985, 1992) and Rudnicki (2002). The analytical results confirmed and extended the findings of earlier, numerical studies.

The aims of the current paper are four-fold: (1) gaining further insight in the poro-mechanical aspects of induced seismicity in displaced faults; (2) clarifying some of the theory underlying (semi-) analytical techniques for fault modelling; (3) investigating the scope for fast (Monte Carlo) simulation of induced seismicity; and (4) developing a semi-analytical framework for embedded fault modelling in large-scale numerical simulation tools. We focus on a semi-analytical approach to describe depletion-induced aseismic fault slip and the onset of seismicity with the aid of Cauchy-type singular integral equations and the use of Chebyshev polynomials for their solution. The core of the underlying theory was developed decades ago by Bilby & Eshelby (1968) and Rice (1968, 1980) on the basis of dislocation theory and fracture mechanics, with the more recent book by Segall (2010) providing an inspiring overview and further development of much of this material. Also, the stability criterion for aseismic slip developed by Uenishi & Rice (2003) forms important input for our paper. We build on this existing body of knowledge to cope with the particularities of the stress field resulting from depletion of a reservoir with displaced faults; these include singularities and discontinuities in the shear and normal stresses, and the development of two distinct slip patches that may merge at increased depletion levels. Our paper also builds on the findings of more recent (numerical) work on depletion-induced seismicity, in particular the seminal reports by Van den Bogert (2015, 2018) and further publications by, amongst others, Buijze et al. (2017, 2019), Buijze (2020) and Van Wees et al. (2017).

The paper is organised as follows. In Section 2, we review expressions for the injection-induced or depletion-induced stress field in a displaced fault. We also introduce an illustrative Example that will be used throughout the paper to highlight various elements of the theory. In Section 3, we discuss aspects of dislocation theory and its relationship to slip-induced stresses, whereafter Section 4 treats the use of Cauchy-type singular integral equations to compute the induced slip pattern for a given depletion-induced stress field. Section 5 illustrates the computation of induced fault slip for the Example corresponding to three different friction models: zero friction, constant static friction and slip-weakening friction. For the latter, we also introduce a modified Uenishi & Rice (2003) stability criterion. Section 6 treats computational aspects, and Section 7 presents a discussion with conclusions about the potential of the developed theory for computational applications. The Appendices provide all mathematics necessary to derive the expressions in the main text. Further mathematical background material is provided in the Supporting Information.

Induced elastic stresses in a displaced fault

Initial and incremental stresses

We consider a two-dimensional (2D) plane-strain model of an infinitely wide reservoir of height $h = a + b$ intersected by a displaced non-sealing normal fault with throw ${t_{f}} = b - a$ ; see Fig. 1 with typical values given in Table 1. We assume the presence of an initial regional stress pattern with principal stresses $\sigma _{yy}^0$ (vertical) and

(1) $$\eqalign{\sigma _{xx}^0 &= {\sigma} _{xx}^{\prime{0}} - \alpha {p^0} \cr &= {K^0}\sigma _{yy}^{\prime{0}} - \alpha {p^0} \cr &= {K^0}(\sigma _{yy}^0 + \alpha {p^0}) - \alpha {p^0},0}$$

Figure 1. Infinitely wide reservoir with a displaced normal fault.

Table 1. Reservoir properties and fault geometry for the Example.

Note: the initial vertical stress, initial pressure and initial effective normal stress have been computed as: $\sigma _{yy}^0(y) = [(1 - \phi ){\rho _s} + \phi {\rho _{fl}}]g(y - {D_0}),\;{\rm{where}}\;\sigma _v^0 \lt\, 0$ , ${p^0}(y) = p_0^0 - {\rho _{fl}}{\kern 1pt} g{\kern 1pt} y$ , $\sigma \prime_ \bot ^0(y) = \sigma _ \bot ^0(y) + \beta {p^0}(y).$ (Valid for reservoir, overburden and underburden).

(horizontal), where $\alpha $ is Biot’s coefficient, ${p^0}$ is the initial pore pressure (a superscript ‘0’ means ‘initial’), ${K^0}$ is the initial effective stress ratio, and where a primed stress variable $\sigma \prime$ represents an ‘effective stress’. We employ a sign convention where positive strains and stresses imply extension and tension, and where pore pressures are positive. The resulting initial normal and shear stresses acting on the fault follow from a coordinate rotation as

(2) $$\sigma _ \bot ^0 = \sigma _{\hat y\hat y}^0 = \sigma _{xx}^0\mathop {\sin }\nolimits^2 \theta + \sigma _{yy}^0\mathop {\cos }\nolimits^2 \theta ,$$

(3) $$\sigma _\parallel ^0 = - \sigma _{\hat x\hat y}^0 = (\sigma _{xx}^0 - \sigma _{yy}^0)\sin \theta \cos \theta ,$$

where $\hat x$ and $\hat y$ are rotated coordinates, and where $\theta $ is the dip angle of the fault; see Fig. 1. A positive-valued shear stress ${\sigma _\parallel }$ corresponds to a normal faulting regime, or in other words, to a situation where the hanging wall (to the left of the fault in Fig. 1) has a tendency to slide down from the foot wall (to the right of the fault). The initial effective normal stress acting at the fault follows as

(4) $$\sigma \prime_ \bot ^0 = \sigma _ \bot ^0 + \beta {p^0},$$

where $\beta $ is an effective stress coefficient which is not necessarily identical to $\alpha $ and is often taken as unity (Scholz, 2019; Fjaer et al., 2021).

An increase or decrease in pore pressure in the reservoir will result in incremental normal and shear stresses in the reservoir and its surroundings. We restrict our analysis to the case of a quasi-steady state, or in other words, to the case of a spatially homogeneous incremental pore pressure $p(t)$ that is a slow function of time t. Closed-form analytical expressions for the corresponding incremental normal and shear stresses in the fault were obtained by Jansen et al. (2019) and can be expressed as

(5) $${\sigma _ \bot } = ({ - {\sigma _{xy}}\sin \theta \cos \theta + {\sigma _{xx}}\mathop {\sin }\nolimits^2 \theta }),$$

(6) $${\sigma _\parallel } = ({{\sigma _{xy}}\mathop {\sin }\nolimits^2 \theta + {\sigma _{xx}}\sin \theta \cos \theta } ),$$

where ${\sigma _{xx}} = {\sigma _{\hat y\hat y}}$ and ${\sigma _{xy}} = - {\sigma _{\hat x\hat y}}$ are normal and shear stresses for a vertical fault, that is, for a dip angle $\theta = {{\pi } \over {2}}$ . They are defined as (see also Appendix A)

(7) $${\sigma _{xx}} = \left\{ {\matrix{ 0 &{{\rm{if}}} & {y \le - b\;\;{\rm{or}}\;\;b \le y} \cr \hfill { - \pi C} \hfill\ &{{\rm{if}}} &\hfill { - b \lt \,y \le - a\;\;{\rm{or}}\;\;a \le y \lt \,b\,,} \cr{ - 2\pi C}\ &{{\rm{if}}} &{ - a \lt \,y \lt \,a} \cr } } \right.$$

and

(8) $${\sigma _{xy}} = \frac{C}{2}\ln \frac{{{{(y - a)}^2}{{(y + a)}^2}}}{{{{(y - b)}^2}{{(y + b)}^2}}},$$

where $C(p(t))$ is a pressure-dependent scaling parameter, with SI units Newton per meter squared, defined as

(9) $$C = {{{(1 - 2\nu )\alpha p(t)}} \over {{2\pi (1 - \nu )}}},$$

with $\nu $ representing Poisson’s ratio, and t time. For dipping as well as vertical faults, the incremental effective normal stress is given by

(10) $${\sigma '_ \bot } = {\sigma _ \bot } + \left\{ {\matrix{ \hfill 0 &{{\rm{if}}}&{y \le - b\;\;{\rm{or}}\;\;b \le y} \cr {\beta p}&{{\rm{if}}}&{ - b \lt \,y \lt\, b} \cr } } \right..$$

In the derivation of Equation 10, it was assumed that only those parts of the fault that are in direct contact with the reservoir experience incremental reservoir pressure, or in other words, that the relevant fault segment is given by $ - b \lt y \lt b$ . If a larger part of the fault is exposed to incremental pressure, the domain where $\beta p$ is added should be extended accordingly.

Relation to other studies

In parallel to the work reported in Jansen et al. (2019), similar analytical expressions were derived by Lehner (2019). Using a slightly different approach, he arrived at a set of alternative closed-form expressions for the induced stresses around a displaced fault, which produce near-identical results when evaluated numerically. More recently, Wu et al. (2021) and Van den Hoek & Poessé (2021) presented similar expressions, based on slightly different derivations. The former publication includes the effect of a (static) pressure difference across the fault, while the latter includes the effect of thermal stresses in addition to incremental pressures.

In all these studies, it is assumed that the reservoir is embedded in an infinite full space, that is, an infinite domain in both horizontal and vertical directions. In reality, the presence of the earth’s surface suggests that the use of an infinite half space might be more appropriate. However, for a small ratio of reservoir thickness over depth, the differences between a full-space and a half-space representation of the stresses are negligible, as demonstrated in detail by Lehner (2019).

At some horizontal distance away from the fault, say at more than a few times the height of the reservoir, the stresses approach those of an infinitely wide reservoir without faults. In the limit $x \to \pm \infty $ , the expressions for an infinitely wide reservoir derived in Jansen et al. (2019) therefore become equal to those for the poro-mechanical behaviour of a straight horizontal ‘elastic thin sheet’ that undergoes unidirectional (vertical) fluid-induced expansion or compression (Bourne & Oates, 2017; Fjaer et al., 2021).

Another, semi-analytical, method to model the elastic stresses in displaced faults was presented by Van Wees et al. (2019) who based their expressions on earlier work by Okada (1992). The underlying solid-mechanical formulation is closely related to those used by Jansen et al. (2019) and Lehner (2019) but the semi-analytical implementation allows, in theory, for the computation of stresses in more complex 3D configurations of fault blocks than the 2D analytical formulations of Jansen et al. (2019), Lehner (2019), Wu et al. (2021) and Van den Hoek & Poessé (2021). Comparison of the stresses in 2D examples using the methods of Van Wees et al. (2019), Lehner (2019) and Jansen et al. (2019) showed near-identical results, while, moreover, all three methods have been checked independently against finite element solutions by their respective authors.

Yet another analytical approach was taken by Hettema (2020) who started from different assumptions leading to different results. However, despite their somewhat different approaches, all the methods referred to in this section share the limitation of assuming elastic material behaviour with uniform homogeneous properties. More importantly, they all disregard the effect of fault slip on the stress field around the fault and the resulting redistribution of stresses. In the present paper, we aim at addressing the latter shortcoming through developing a (semi-) analytical approach that includes the effect of fault slip for various friction laws.

Fault slip and coulomb stress

Fault slip is defined as

(11) $$\delta (s,t) = u_\parallel ^ + (s,t) - u_\parallel ^ - (s,t),$$

where $u_\parallel ^ - $ and $u_\parallel ^ + $ are the along-fault displacements at both sides of the fault; it is governed by combined (i.e. initial plus incremental) shear and effective normal stresses

(12) $${\textstyle \sum}{_\parallel } = \sigma _\parallel ^0 + {\sigma _\parallel },$$

(13) $${\textstyle \sum} {^\prime_ \bot } = \sigma _ \bot ^{\prime{0}} + {\sigma _ \bot^{\prime} }.$$

Slip-provoking conditions occur when

(14) $$|{\textstyle \sum} {_\parallel }| \gt {\textstyle \sum} {_{sl}},$$

where ${{\textstyle \sum} _{sl}}$ is the slip threshold, defined as

(15) $${\textstyle \sum} {_{sl}} = \kappa - \mu {\textstyle \sum}{^\prime_ \bot },$$

with $\kappa \ge 0$ indicating cohesion and $\mu $ the friction coefficient, and where it should be kept in mind that negative normal stresses correspond to compression. In the most general case, the friction coefficient ${\mu ({{\delta}}, {\dot{\delta}} ,\psi ,s,t)}$ is lithology-dependent and thus a function of the along-fault coordinate s. Moreover, it is a function of the slip rate

(16) $${{\dot{\delta}}} (s,t) ={{{d\delta (s,t)}} \over {{dt}}};$$

the cumulated slip

(17) $${{{\tilde{\delta}}}} = \int_0^t |{{\dot{\delta}}} (s,t)|dt;$$

possibly one or more additional state variables ${\psi _i}(s,t)$ and time t directly. Equation 14 implies that slip of the hanging wall may occur in upward or downward direction, where exceedance of the slip threshold ${{\textstyle \sum} _{sl}}$ by a positive combined shear stress ${{\textstyle \sum} _\parallel }$ implies downward slip of the hanging wall, or in other words, a continued normal fault development. In the remainder of the paper, we will only consider such downward slip without reversals of direction. Therefore, we have $\,\,{\tilde{\delta} = \delta}$ , and we can employ the usual definition of the pre-slip Coulomb stress

(18) $${\textstyle \sum} {_C}(s,t) = {\textstyle \sum} {_\parallel }(s,t) - {\textstyle \sum} {_{sl}}(s,t),$$

in which slip corresponds to positive values of ${{\textstyle \sum} _C}$ .

Illustrative example

Figure 2 (left) depicts an example of combined shear stresses and the (downward) slip threshold for a displaced fault crossing a hydrocarbon reservoir that is experiencing a gradual pore pressure decrease due to quasi-steady-state depletion. The parameter values for this example have the same order of magnitude as those of the Groningen natural gas reservoir in the Netherlands (NAM, 2016), and their values are listed in Table 1 together with calculation details of the initial stresses and pore pressures. We will use this Example throughout the paper to illustrate various elements of the theory.

Figure 2. Shear stresses, slip threshold and pre-slip Coulomb stresses for the Example. Left: red lines represent the combined depletion-induced shear stresses ${\textstyle \sum} {_\parallel }$ and black lines the slip threshold ${\textstyle \sum} {_{sl}}$ as a function of vertical position $y$ . Right: red lines represent the pre-slip Coulomb stresses ${\textstyle \sum} {_C} = {\textstyle \sum} {_\parallel } - {\textstyle \sum} {_{sl}}$ . In both the left and right figures the green areas indicate those regions where the shear stresses exceed the slip threshold and, accordingly, the green horizontal lines at $y = - 76,{\kern 1pt} - 51,{\kern 1pt} 47$ and $76$ m correspond to the intersections between shear stresses and slip threshold, i.e. the zeros of the pre-slip Coulomb stresses. The light gray areas depict the vertical positions of the foot wall and the hanging wall, and the dash-dotted rectangle in the right figure corresponds to the detailed view in Figure 5 (left).

The figure corresponds to an incremental pressure $p = - 25$ MPa (where the negative sign implies depletion) and displays two potential slip patches formed by areas where the shear stress exceeds the slip threshold, indicated in green. The two patches are initiated at the ‘internal’ reservoir-fault corners, at $y = \pm a = \pm 75$ m, and, with increasing depletion, gradually grow into the reservoir such that the top patch extends downwards and the bottom patch upwards. Figure 2 (right) displays the corresponding Coulomb stresses with green areas for the same y values as in the left figure. Note that these results do not yet include the effects of slip.

With continuing depletion, the two potential slip patches may merge and form a single patch. This stress pattern development is typical for depletion-induced shear stresses in displaced faults and was first described by Van den Bogert (2015, 2018) and subsequently by others; see Buijze et al. (2017, 2019), Van Wees et al. (2017, 2018, 2019) for numerical studies, and Jansen et al. (2019) and Lehner (2019) for analytical approaches. Opposedly, injection-induced shear stresses result in the development of potential slip patches at the ‘external corners’ (in this example located at $y = \pm 150$ m) which subsequently (mainly) grow outwards into the overburden and underburden; see Jansen et al. (2019).

Intersections

Intersections of the shear stress profiles with the slip threshold in Fig. 2 (left) have been indicated with horizontal green lines. They are identical to the zeros of the Coulomb stresses; see Fig. 2 (right). Their magnitudes can be obtained by solving iteratively for y from the implicit equation

(19) $${\textstyle \sum} {_\parallel }(y) = {\textstyle \sum} {_{sl}}(y),$$

with the aid of Equations 1–15 resulting in four values $\{ {y_1},{y_2},{y_3},{y_4}\} $ , where

(20) $$ - b \lt \,{y_1} \lt - \!a \lt \,{y_2} \lt \,0 \lt \,{y_3} \lt \,a \lt \,{y_4} \lt \,b,$$

as long as the slip patches have not merged, and two relevant values $\{ {y_1},{y_4}\} $ , where

(21) $$ - b \lt \,{y_1} \lt - \!a \lt \,0 \lt \,a \lt \,{y_4} \lt \,b,$$

thereafter. Note: in Fig. 2, there are four intersections, but the positions of ${y_1}$ and ${y_4}$ are just below $y = - a$ and just above $y = a$ , respectively, such that they nearly coincide with the dotted black lines that indicate the locations of the top of the hanging wall and the bottom of the foot wall.

We emphasise that the stress distribution as depicted in Fig. 2 is based on an elastic equilibrium without taking into account the effect of slip. Once slip occurs, whether aseismic or seismic, the distribution will change, as will be shown below, and so will the intersections. These pre-slip intersections should therefore not be used as a basis for a quantitative prediction of slip patch size, but rather as an indication of the potential region where slip may occur.

Logarithmic singularities and jump discontinuities

Figure 2 (left) displays sharp peaks in the shear stresses and the slip threshold at $y = \pm a$ and $y = \pm b$ , and Fig. 2 (right) shows similar peaks in the pre-slip Coulomb stresses. These peaks correspond to singularities in the underlying logarithmic expressions which result in infinite values when evaluated for $y = \pm a$ or $y = \pm b$ ; see Equations 5 and 6 which both inherit the logarithmic term from Equation 8 for the shear stress in a vertical fault. Equations 5 and 6 also both contain jump discontinuities at $y = \pm a$ and $y = \pm b$ . Although these are hidden in Fig. 2 by the peaks, their presence can be inferred from their common descendance from the expression for the normal stresses in a vertical fault; see Equation 7 which displays jump discontinuities of magnitude $C\pi $ at $y = \pm a$ and $y = \pm b$ .

A consequence of the logarithmic singularities is that numerical results, for example from finite element simulations, that suggest the existence of a ‘threshold depletion’ above which there seems to be an absence of slip, are mesh-size dependent and may fail to converge to a fixed value on refinement of the simulation grid. In our model, the logarithmic singularities as well as the jump discontinuities result from the sharp ‘internal’ and ‘external’ reservoir-fault corners at $y = \pm a$ and $y = \pm b$ . In reality, such infinite peaks are physically impossible, and small amounts of aseismic fault slip, plastic deformation and/or pore pressure diffusion will result in finite peak stresses, while the presence of more rounded corners or gradual pore pressure transitions between the reservoir and the surrounding rock will have the same result. Mathematically, such a smoothing of the peaks can be obtained through regularising the expressions for the stresses; see Appendix B. However, even under such more realistic conditions, a major effect of the presence of a displaced fault in a depleting reservoir is the occurrence of high-magnitude stress peaks.

Dislocation theory and singular integral equations

Dislocations and double couples

To obtain insight into the relation between depletion-induced stresses and fault slip, we use elements of dislocation theory as developed in the field of material sciences (Burgers, 1939; Nabarro, 1952; Weertman, 1996; Cai & Nix, 2016). Since the 1950s, dislocations have been applied to represent earthquake sources by, for example, Vvedenskaya (1956), Keylis-Borok (1956) and Steketee (1958a, b). We refer to Eshelby (1973), Segall (2010), Udías et al. (2014) and references therein for further discussions on the use of dislocation theory for geophysical applications.

An edge dislocation (also known as a glide dislocation) can be represented in $(x,y)$ coordinates as an in-plane shear displacement along a semi-infinite slip line; see Fig. 3 (top left) where the material just to the right of the y axis in the half plane $\{ y \lt 0\} $ has been displaced in the positive y direction over a distance ${\overline{u}_y} = \delta /2$ while the material in the same half plane but just to the left of the y axis has been displaced in the negative y direction over the same distance. As a result, the displacement field contains a singularity at the origin.

Figure 3. Top left: slip $\delta $ (red) and vertical displacements ${\overline{u}_y}$ (orange) in a vertical fault corresponding to an edge dislocation at the origin with given displacements $\overline{u}_y^ - $ and $\overline{u}_y^ + $ along the half line $\{ x = 0,y \lt 0\} $ . The corresponding slip $\delta = \overline{u}_y^ + - \overline{u}_y^ - $ has a magnitude of 0.05 m. Top right: slip-induced shear stresses ${\overline \sigma _\parallel }$ in the fault (which are identical to ${\overline \sigma _{xy}}$ in this vertical case) resulting from the dislocation. Bottom left: slip $\delta $ (red) and vertical displacements ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over u} _y}$ in a vertical fault corresponding to a ramped slip profile between boundaries ${y_ - } = - 100$ m and ${y_ + } = 100$ m with slope $c = -{{{0.05}} \over {{200}}}$ . Bottom right: slip-induced shear stresses ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel } = {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _{xy}}$ resulting from this ramped profile. The material properties to compute the stresses in the sub figures at the right have been chosen according to Table 1.

Expressions for the plane-strain displacements ${\overline{u}_i}$ and stresses ${\overline \sigma _{ij}},i \in \{ x,y\} ,j \in \{ x,y\} $ , around an edge dislocation, can be obtained with techniques from the theory of elasticity (Nabarro, 1952; Barber, 2010; Cai & Nix, 2016). Alternatively, one can use the fact that an edge dislocation with slip of magnitude $\delta = \overline{u}_y^ + - \overline{u}_y^ - $ results in a seismic moment per unit area of magnitude $m = \delta G$ , where G is the shear modulus. In 2D, the moment at a given point can then be interpreted as a nucleus of strain in the form of a double couple $\{ m, - m\} $ per unit length where the positive and negative couples act counter-clockwise and clockwise, respectively, a well-known concept in earthquake source mechanics. For further information, see, for example, Steketee (1958b), Eshelby (1973), Segall (2010) and Udías et al. (2014). Section S1 of the Supporting Information describes how to obtain expressions for ${\overline{u}_i}$ and ${\overline \sigma _{ij}}$ , using Green’s functions for a point source (Lord Kelvin, 1848) and the nucleus-of-strain concept (Love, 1927) to represent double couples. For an inclined fault passing through the origin, and with an edge dislocation just there, it follows that the stresses at the fault location can be expressed as

(22) $${\overline \sigma _ \bot } = 0,$$

(23) $${\overline \sigma _\parallel } = {{{\delta G}} \over {{2\pi (1 - \nu )s}}} = {{{\delta G\sin \theta }} \over {{2\pi (1 - \nu )y}}},$$

where s is the along-fault coordinate, and $\delta $ is the along-fault slip for $s \lt 0$ (see Fig. 1). Figure 3 (top right) displays the slip-induced shear stresses ${\overline \sigma _\parallel }$ for the case of an edge dislocation in a vertical fault in which case it holds that ${\overline \sigma _\parallel } = {\overline \sigma _{xy}}$ . Note that Equation 22 implies that along-fault displacements caused by slip do not result in a change in normal stresses in that fault and therefore do not directly influence the magnitude of the slip thresholds.

Distributed dislocations

Next we consider an edge dislocation with an along-fault slip $\delta $ that is no longer constant but is a function of the along-fault coordinate $s(y) = {{y} \over {{\sin \theta }}}$ . The shear stresses ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }$ resulting from an array of infinitesimal edge dislocations $d\delta $ with magnitude

(24) $$d\delta = {{{\partial \delta (s)}} \over {{\partial s}}}ds = \nabla \delta (s)ds,$$

for ${s_ - } \le s \le {s_ + }$ and $d\delta = 0$ otherwise, can be expressed with the aid of Equation 23 as

(25) $${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }(s) = A\int_{ - \infty }^\infty {{{\nabla \delta (\xi )}} \over {{\xi - s}}}\,d\xi = A\int_{{s_ - }}^{{s_ + }} {{{\nabla \delta (\xi )}} \over {{\xi - s}}}\,d\xi ,$$

where, for plane-strain conditions,

(26) $$A = {{G} \over {{2\pi (1 - \nu )}}},$$

and

(27) $$\nabla \delta (\xi ) = {\left. {{{\partial \delta (s)}} \over {{\partial s}}} \right|_{s = \xi }};$$

see also Bilby & Eshelby (1968). The quantity $\nabla \delta (s)$ can be interpreted as the dislocation density or the slip gradient.

The integral in Equation 25 is a so-called Cauchy-type singular integral of the first kind or a Cauchy integral for short (Cates, 2019). This means that, although it contains a singularity at $\xi = s$ , it is associated with a finite result that is known as its principal value (PV) defined as

(28) $${\rm{PV}}\int_{{s_ - }}^{{s_ + }} {{{\nabla \delta (\xi )}} \over {{\xi - s}}}\,d\xi = \mathop {\lim }\limits_{\varepsilon \downarrow 0} \left\{ {\int_{{s_ - }}^{s - \varepsilon } {{{\nabla \delta (\xi )}} \over {{\xi - s}}}\,d\xi + \int_{s + \varepsilon }^{{s_ + }} {{{\nabla \delta (\xi )}} \over {{\xi - s}}}\,d\xi } \right\}.$$

These integrals frequently occur in contact and fracture mechanics and a classic mathematical text that extensively treats theoretical and applied aspects is the book by Muskhelishvili (1953). Somewhat more accessible treatments are given by, for example, Tricomi (1957), Weertman (1996) and Estrada & Kanwal (2000). The first integral in Equation 25 can also be interpreted as an (infinite) Hilbert transform as used in signal processing. Applications in geophysics were pioneered in the 1960s by Bilby & Eshelby (1968) and Rice (1968), with a more recent treatment being given in the book by Segall (2010). In the remainder of the paper, we will not explicitly indicate the PV for singular integrals.

If the slip gradient $\nabla \delta (s)$ is a simple function of s, it may be possible to obtain the PV through straightforward analytical integration. As an example, the step-shaped slip profile in the vertical fault in Fig. 3 (top left) can be replaced by a ramped profile by choosing a constant slip gradient $\nabla \delta (y) = c \lt 0$ over a slip interval ${y_ - } \le y \le {y_ + }$ , where ${y_ \pm } = \pm 100$ m, and a zero gradient otherwise. (Note that in a vertical fault we have $s = y$ ). A possible realisation of the corresponding vertical displacements $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over u} _y^ - $ and $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over u} _y^ + $ has been depicted in Fig. 3 (bottom left), and the resulting shear stress follows from Equation 25 as

(29) $${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel } = {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _{xy}} = {{{Ac}} \over {2}}\ln \left[ {{{{{(y - {y_ + })}^2}}} \over {{{{(y - {y_ - })}^2}}}} \right],$$

and is depicted in Fig. 3 (bottom right). Note that it is only the slip gradient $\nabla \delta $ that determines the slip-induced stresses and not the magnitude $\delta $ of the slip itself. Therefore, an arbitrary constant may be added to $\delta $ , over the entire infinite domain, and the same holds for the corresponding displacements.

It can be seen that the peaks in the stress profile at $y = {y_ \pm } = \pm 100$ m in Fig. 3 (bottom right) resemble those in the shear stress profile in Fig. 2 (left). This resemblance can also be noted in the logarithmic terms in Equations 8 and 29, a feature that will be made use of later on in this paper.

Computing the slip gradient and the slip

Mathematical aspects

Equation 25 provides the shear stress ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }(s)$ for a given slip gradient $\nabla \delta (s)$ . However, in the following we will consider situations where we aim to anihilate the pre-slip Coulomb stress ${{\textstyle \sum} _C}(s)$ over the slip patches through finding a slip distribution $\delta (s)$ such that ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }(s) = - {{\textstyle \sum} _C}(s)$ . In that case, we need to compute the slip gradient for a given shear stress which requires the inverse of Equation 25. It was shown in Muskhelishvili (1953) that such an inverse can be obtained analytically if the function ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }(\xi )$ is Hölder continuous, where Hölder continuity is a stricter form of continuity than the version used in regular mathematical analysis; see the Supporting Information belonging to this paper. As discussed before, Equations 5 and 6 for the shear and normal stresses in a displaced fault display jump discontinuities in $y = \pm a$ and $y = \pm b$ . (An exception is a vertical fault in which case the shear stresses are continuous but the normal stresses still display jump discontinuities). Therefore, these stresses are not even continuous in the regular sense in those points, leave alone Hölder continuous, and the same holds for the pre-slip Coulomb stresses ${{\textstyle \sum} _C}$ .

A more detailed analysis of the discontinuities reveals that inversion results in slip values that contain a weak (logarithmic) singularity in the expression for the slip gradient $\nabla \delta $ , while it produces a (just) regular result for the slip $\delta $ ; see the Supporting Information. However, as discussed before, we can regularise Equations 5 and 6 such that the jump discontinuities and logarithmic singularities are removed; see also Appendix B. In that case, we can safely use the approach of Muskhelishvili (1953) to invert Equation 25 which results in a solution that is well-known in the fracture mechanics and contact mechanics literature as will be discussed briefly below and in more detail in Appendix C.

Inverse equation

For the case that $\nabla \delta $ remains finite at both $s = {s_ - }$ and $s = {s_ + }$ , it was shown by Muskhelishvili (1953) that the inverse of Equation 25 is given by

(30) $${\nabla \delta (s) = {{{\Psi (s)}} \over {{{\pi ^2}A}}}\int_{{s_ - }}^{{s_ + }} {{{{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } }_\parallel }(\xi )}} \over {{\Psi (\xi )\,(s - \xi )}}}\,d\xi ,\;{s_ - } \le s \le {s_ + },}$$

where

(31) $$\Psi (s) = \sqrt {(s - {s_ - })({s_ + } - s)} ,$$

with a similar expression for $\Psi (\xi )$ . Equation 30 is only valid if the following two conditions are fulfilled (Bilby & Eshelby, 1968):

(32) $$\int_{{s_ - }}^{{s_ + }} {{{\xi ^k}\,{{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } }_\parallel }(\xi )}} \over {{\Psi (\xi )}}}\,d\xi = 0,\;k = 0,1.$$

For various proofs, see Muskhelishvili (1953), Weertman (1996), Segall (2010) and the Supporting Information belonging to this paper. For a given function ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }$ , conditions (32) can only be met for specific values of ${s_ - }$ and ${s_ + }$ . These two values should therefore be solved for from conditions (32), which typically requires an iterative process. Physically, conditions (32) imply that the sum of the slip-induced shear stresses as well as their first moment vanish when evaluated over the slip patch ${s_ - } \lt \,s \lt \,{s_ + }$ . In terms of the corresponding slip gradient, this means that the magnitudes of the gradient at the endpoints are finite, in other words, that the slip at those points is equal to zero.

Once the slip gradient $\nabla \delta (s)$ has been established with Equation 30, the along-fault slip can be determined with the aid of Equation 27 as

(33) $$\delta (s) = \int_{{s_ - }}^s \nabla \delta (\xi )\,d\xi ,$$

where we can use our knowledge that $\delta ({s_ - }) = 0$ . It should be noted that Equation 25 is valid for $ - \infty \lt \,s \lt \,\infty $ , whereas the validity of Equation 30 is restricted to ${s_ - } \le s \le {s_ + }$ such that the latter is now equivalent to a finite Hilbert transform (Tricomi, 1957; Weertman, 1996).

Equations 30 and 32, with ${s_ - } = - {s_ + }$ , were used in Bilby & Eshelby (1968) and Rice (1968) to compute the stresses in various fracture configurations modelled as arrays of dislocations, and in Rice (1980) and Segall (2010) in a similar fashion to compute stresses in and around faults caused by imposed shear stresses. In Uenishi & Rice (2003), the equations served as a basis to determine the stability of a fault loaded by a peaked shear stress. In the following, we will start from Equations 30 and 32 to obtain the full slip distribution along a displaced fault and thereafter return to the stability aspects.

Slip-induced shear stresses and slip gradient

From now on, we will express all variables, including the along-fault slip $\delta $ and the along-fault slip gradient $\nabla \delta = {{{d\delta }} \over {{ds}}}$ as functions of the vertical coordinate $y = s\sin \theta $ . As discussed before, a decrease in pore pressure in a reservoir with a displaced fault will always result in the development of two potential slip patches in those fault segments where the shear stresses ${\textstyle \sum} {_\parallel }(y,t)$ exceed the slip threshold ${\textstyle \sum} {_{sl}}(y,t)$ , that is, where the pre-slip Coulomb stresses ${\textstyle \sum} {_C}(y,t)$ are positive; see Fig. 2. A complicating factor is that the slip in one patch influences the shear stresses in the other patch and vice versa, although we may disregard this coupling effect as long as the patches are located far enough from each other.

In any case, we seek a slip gradient $\nabla \delta (y,t)$ , and thus a slip pattern $\delta (y,t)$ , that results in a slip-induced shear stress distribution ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }(y,t)$ that annihilates the pre-slip Coulomb stress defined in Equation 18, such that the green areas in Fig. 2 vanish. Recall that slip does not influence the normal stresses; see Equation 22. However, such a slip pattern also induces a change in shear stresses for y values outside those green fault segments and therefore influences the location of the intersections ${y_i}$ . As a result, the region where the slip-induced stresses annihilate the pre-slip Coulomb stresses becomes larger. Thus, as long as the coupling effect can be disregarded, we are looking for a slip pattern with corresponding slip patch boundaries that obey the following mixed boundary conditions for ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }$ and $\nabla \delta $ :

(34) $${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }(y,t) = - {\textstyle \sum} {_C}(y,t),\,\,\,\{\,\tilde{y}_1(t) \lt y \lt \tilde{y}_2(t)\} \vee \{\tilde{y_3}(t) \lt y \lt \tilde{y_4}(t)\},$$

(35) $$\nabla \delta (y,t) = 0,\;\;\{ \,y \le\, {{{\tilde{y}}}_1}(t)\} \vee \{ {{{\tilde{y}}}_2}(t) \le y \le {{{\tilde{y}}}_3}(t)\} \vee \{ {{{\tilde{y}}}_4}(t) \le y\} ,$$

where $\,{\tilde{y}_i}(t),i = 1, \ldots ,4$ are shifted intersections that now constitute the post-slip patch boundaries, and where ${\textstyle \sum} {_C}(y,t)$ is given by Equation 18 with further details in Equations 1–15. Here we explicitly indicate the time dependency of the variables which may result from one or more sources, such as the slow (i.e. quasi-steady-state) change in reservoir pressure $p(t)$ , the time dependency of the cohesion $\kappa (y,t)$ or the rate and state dependency of the friction coefficient ${\mu (y,t) = \mu (\tilde{\delta} (y,t),\ {\dot \delta} (y,t),{\psi _i}(y,t))}$ . From now on, we will refrain from indicating these dependencies except when necessary. It is shown in Appendix C how Equations 30–32 can be combined with Equations 34 and 35 to obtain a system of equations to solve for the boundaries $\{ {y_ - },{y_ + }\} $ and a Cauchy integral formulation to compute the slip gradient $\nabla \delta (y)$ in each of the two slip patches.

For increasing depletion, both slip patches will grow in size and approach each other in which case the coupling effect may no longer be disregarded. In that case, we can still use Equations 30–32 in combination with 34 and 35 if we replace the term $ - {\textstyle \sum} {_C}(y)$ in Equation 34 by $ - {\textstyle \sum} {_C}(y) - \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel ^ \times (y)$ , where ${\textstyle \sum} {_C}(y)$ are now the pre-slip Coulomb stresses in a patch while the cross-term $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel ^ \times (y)$ represents the shear stresses in that same patch as caused by slip in the other patch; see Appendix C for further details. Alternatively, we could employ a split-integral approach, using results from fracture mechanics for multiple dislocation zones, resulting in two two-term integrals to compute the slip gradients in the two slip patches (Weertman, 1996). The integrals then contain a square root of a fourth-order polynomial, resulting in a set of four conditions that should be fulfilled to obtain the four slip patch boundaries. Although this alternative formulation has the benefit of avoiding iterations, it requires a simultaneous solution of four highly nonlinear equations. For the examples that we considered, this turned out to be much more prone to generating spurious results than the iterative sequential solution of two sets of two nonlinear equations, and we therefore did not pursue the split-integral approach any further.

After merging of the slip patches, which may occur at larger depletion values, the coupling effect is no longer relevant. In that case, the slip and the slip patch boundaries can be computed with the aid of the same equations as used for the uncoupled patches, but now applied to a single merged patch.

Numerical and semi-analytical integration

Evaluation of the Cauchy integrals required to compute the slip gradient can be performed numerically if attention is paid to handling of the singularities; see, for example, Golberg (1990), Keller & Wróbel (2016) and Viesca & Garagash (2018) who provide several algorithms that benefit in varying degrees from the underlying mathematical structure of the equations. Alternatively, semi-analytical expressions can be obtained by expanding the known function in the integrand in terms of Chebyshev polynomials (Mason & Handscomb, 2003; Uenishi & Rice, 2003; Segall, 2010; Barber, 2018). We applied such a semi-analytical approach for the various Cauchy integrals, with details as discussed in Appendix D. To verify the semi-analytical results, we also obtained results through numerical integration with a relatively simple ‘staggered grid’ approach inspired by the paper of Uenishi & Rice (2003). Further details are provided in Appendix E, which also contains a brief comparison of the computational performance of the two approaches.

The influence of friction

No friction

As a first step in addressing the effects of friction on the development of slip, consider the hypothetical case of a frictionless fault. This implies a restriction to vertical faults because the absence of friction would make it impossible to sustain the initial shear stresses corresponding to the usual situation with unequal horizontal and vertical principal stresses. For this case of a vertical frictionless fault without initial shear stresses, Bourne & Oates (2017) developed an approximate solution in which they assumed that, in each of the two fault blocks, the post-slip displacements become uniaxial (vertical) with values equal to those at $x = \pm \infty $ . The same assumption, following a slightly different derivation, was made in the Supporting Information of Jansen et al. (2019) leading to the same approximate solution. Here we relax this assumption and start from the stresses in a vertical fault which were given in Equations 7 and 8.

For such a vertical frictionless fault, we do not yet need an inverse formulation but can directly apply Cauchy Equation 25 for the Coulomb stresses ${\textstyle \sum} {_C}$ as function of a given slip gradient $\nabla \delta $ as follows. We know that the pre-slip Coulomb stresses are just equal to the incremental shear stresses and are therefore given by Equation 8. Now recall the resemblance between the stresses induced by depletion and those resulting from a ramped slip gradient, see Equations 8 and 29. Exploiting this analogy, and guided by the approximate solution of Bourne & Oates (2017), we arrive at the following expression for the depletion-induced slip in a frictionless vertical fault:

(36) $$\delta (y) = {{\gamma C} \over A} \cdot \left\{ {\matrix{ 0 &{{\rm{if}}}&&y&{ \le\! - b,} \cr { - (y + b)}&{{\rm{if}}}& { - b\! \lt }& y &{ \le\! - a,} \cr {\;\;{\kern 1pt} (a - b)}& {{\rm{if}}}&{ - a\! \lt } &y &{ \lt \,a,} \cr {\;\;{\kern 1pt} (y - b)}& {{\rm{if}}}& {{a} \!\le } &y &\,{ \lt b,} \cr 0 &{{\rm{if}}}& {{b} \le } &y &, \cr } } \right.$$

where $\gamma $ is an auxiliary variable which is equal to one for our proposed solution while it is equal to ${{1} \over {{2(1 - \nu )}}}$ for the approximate solution of Bourne & Oates (2017). The shear stresses in the fault corresponding to Equation 36 can be determined with the aid of Equation 25 as

(37) \begin{align} {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }(y) & = \gamma C\left( {\int_{ - b}^{ - a} {{{ - 1}} \over {{\xi - y}}}\,d\xi + \int_a^b {{1} \over {{\xi - y}}}\,d\xi } \right) \nonumber \\ & = - {{{\gamma C}} \over {2}}\ln \left[ { {{{{{(y - a)}^2}{{(y + a)}^2}}} \over {{{{(y - b)}^2}{{(y + b)}^2}}}}} \right].\end{align}

For $\gamma = 1$ , these slip-induced shear stresses are exactly equal to minus the depletion-induced shear stresses ${\sigma _\parallel }$ given in Equation 8 such that we find for the post-slip Coulomb stresses:

(38) $${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over {\textstyle \sum} } _C} = {\textstyle \sum} {_C} + {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel } = {\sigma _\parallel } - {\sigma _\parallel } = 0,$$

which proves that Equation 36 is indeed a correct solution for slip in a frictionless vertical fault. We note that other, equally valid solutions could be obtained by adding an arbitrary constant amount of slip ${\delta _0}$ to the slip distribution of Equation 36 over the entire fault range $ - \infty \lt y \lt \infty $ . However, even in case of a infinitesimal amount of friction, the slip at $y = \mp \infty $ would vanish, because the slip-induced shear stresses as given in Equation 8 vanish at $y = \mp \infty $ , and therefore, it must hold that ${\delta _0} = 0$ .

Figure 4 (left) depicts the pre-slip Coulomb stresses ${{\textstyle \sum} _C}$ , the slip-induced shear stresses ${\mathord{\buildrel{\lower-5.5pt\hbox{$\!\!\!\!\!\scriptscriptstyle\smile$}} \over {\textstyle \sum} } _\parallel}$ and the resulting post-slip Coulomb stresses ${\mathord{\buildrel{\lower-5.5pt\hbox{$\!\!\!\!\!\scriptscriptstyle\smile$}} \over {\textstyle \sum} } _C} = {{\textstyle \sum} _C} + {\mathord{\buildrel{\lower-5.5pt\hbox{$\!\!\!\!\!\scriptscriptstyle\smile$}} \over {\textstyle \sum} } _\parallel }$ for a modified version of the Example, with $\theta = 90$ deg. and $\mu = 0$ such that it now represents a frictionless vertical fault. The solid red line in Fig. 4 (right) depicts the corresponding fault slip $\delta $ as given in Equation 36 with $\gamma = 1$ . The dashed line represents the approximate solution of Bourne & Oates (2017) which is given by Equation 36 with $\gamma = {{1} \over {{2(1 - \nu )}}}$ . We verified our semi-analytical solution against a fully numerical one obtained with an in-house finite-volume code for coupled porous media flow and mechanics (Novikov et al., 2022). The numerical results are indicated with blue triangles in Fig. 4 (right) and they almost exactly coincide with the solid red line.

Figure 4. Depletion-induced stresses and fault slip in a frictionless vertical fault (modified Example). Left: pre-slip Coulomb stresses ${{\textstyle \sum} _C}$ (red), slip-induced shear stresses ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }$ (orange) and post-slip Coulomb stresses ${\mathord{\buildrel{\lower-5.5pt\hbox{$\!\!\!\!\!\scriptscriptstyle\smile$}} \over {\textstyle \sum} } _C}$ (blue dots). The green area indicates the region where the shear stresses exceed the (zero) slip threshold and, accordingly, the green horizontal lines at $y = \pm 119$ m correspond to the zeros ${y_1}$ and ${y_4}$ of the pre-slip Coulomb stresses. Right: the solid red line represents the corresponding fault slip $\delta $ . The dashed red line represents the approximate solution of Bourne and Oates (2017). The blue triangles depict numerical results obtained with a finite volume code. The blue horizontal lines at $y = \pm 150$ m correspond to the boundaries ${\tilde{y}_1}$ and ${\tilde{y}_4}$ of the slip patch.

In this section, it was shown that the approximate expression for slip in a frictionless fault as first derived by Bourne & Oates (2017) underestimates the present result with a factor ${\gamma ^{ - 1}} = 2(1 - \nu )$ . The same holds for the corresponding maximum possible seismic moment per unit strike length ${\hat M_s}$ which is now given by

(39) \begin{align}{\hat M_s} &= G\int_{ - b}^b \delta (y)\,dy = - \pi {\gamma ^{ - 1}}C({b^2} - {a^2}) = - \pi {\gamma ^{ - 1}}Ch{t_{f}} \\ &= (2\nu - 1)\alpha ph{t_{f}},\end{align}

where $h = a + b$ and ${t_{f}} = b - a$ are the reservoir height and fault throw. This expression forms an upper boundary to the seismic moment per unit strike length that can be generated by depletion-induced completely seismogenic fault slip. However, it excludes the potential effects of propagation of the fault slip above or below the reservoir.

Constant friction

As a next step in addressing the effects of friction, consider the depletion-induced elastic stresses in the Example with a homogeneous and constant static friction coefficient $\mu (y,t) = {\mu _{st}} = 0.52$ . We now employ the inverse formulation discussed above and given in detail in Appendix C and perform semi-analytical integration using Chebyshev polynomials as described in Appendix D, verified with numerical integration using a ‘staggered grid’ approach as described in Appendix E.

Consider a gradually increasing depletion, in other words, a gradually decreasing incremental pressure $p(t) \lt \ 0$ . If no slip would occur, the pre-slip Coulomb stresses would gradually grow and so would the (potential) slip patches with boundaries ${y_i},i = 1, \ldots ,4$ ; see Equation 20. However, from the results obtained by Uenishi & Rice (2003) we know that, for a constant friction coefficient, stable non-seismic slip will occur. Figure 5 (left) displays the pre-slip Coulomb stresses ${{\textstyle \sum} _C}$ , the slip-induced shear stresses ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }$ (originating from slip in both patches) and the resulting post-slip Coulomb stresses ${\mathord{\buildrel{\lower-5.5pt\hbox{$\!\!\!\!\!\scriptscriptstyle\smile$}} \over {\textstyle \sum} } _C} = {{\textstyle \sum} _C} + {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }$ for the original Example, that is, just like in Fig. 4 but now with $\theta = 70$ deg. and ${\mu _{st}} = 0.52$ . The figure illustrates that, inside the green areas, the negative-valued slip-induced stresses (orange lines) exactly annihilate the positive-valued pre-slip Coulomb stresses (red lines) such that the post-slip Coulomb stresses (blue dots) become just zero. However, slip also results in significant positive-valued shear stresses just above and below the green areas, and as a result, the slip patches will grow. In case of stable non-seismic slip, as is the case for a static friction coefficient without slip weakening or velocity weakening, the growth process will stabilise such that slip patches, that is, regions of zero post-slip Coulomb stresses, will form between boundaries ${\tilde \!\!\!\!\,y_1 \lt y \lt \tilde \!\!\!\!\,y_2}$ and $\tilde \!\!\!\!\,y_3 \lt y \lt \tilde \!\!\!\!\, y_4$ (Uenishi & Rice, 2003). This growth process is also apparent from the green and blue horizontal lines in Fig. 5: the green lines indicate the boundaries ${y_i},i = 1, \ldots 4$ , of the original potential slip patches while the blue lines indicate the boundaries ${\tilde \!\!\!\!\,y_i}$ of the patches once slip has occurred. Figure 5 (right) displays the corresponding fault slip for increasing depletion and illustrates the gradual growth and subsequent merging of the patches.

Figure 5. Depletion-induced stresses and fault slip for the Example. Left: pre-slip Coulomb stresses ${{\textstyle \sum} _C}$ (red), slip-induced shear stresses ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }$ (orange) and post-slip Coulomb stresses ${\mathord{\buildrel{\lower-5.5pt\hbox{$\!\!\!\!\!\scriptscriptstyle\smile$}} \over {\textstyle \sum} } _C}$ (blue dots) for $p = - 25$ MPa. The green areas indicate the regions where the pre-slip shear stresses exceed the slip threshold. Accordingly, the green horizontal lines at $y = - 76,{\kern 1pt} - 51,{\kern 1pt} 47$ and $76$ m (partly hidden by the red curves) correspond to the zeros ${y_i},\,i = 1, \ldots ,4,$ of the pre-slip Coulomb stresses (see also Figure 2). The blue horizontal lines at $y = - 80,{\kern 1pt} - 33,{\kern 1pt} 29$ and $80$ m correspond to the boundaries ${\tilde \!\!\!\!\,y_i}$ of the slip patches. The black curves represent the pre-slip Coulomb stresses ${{\textstyle \sum} _C}$ for other depletion values ( $p = - 24,{\kern 1pt} - 26,{\kern 1pt} - 27$ and $ - 28$ MPa). Right: Corresponding fault slip showing the transition from two slip patches to a single, merged slip patch for increasing depletion. The red curve corresponds to a depletion pressure of $ - 25$ MPa, and the black curves to $p = - 24,{\kern 1pt} - 26,{\kern 1pt} - 27$ and $ - 28$ MPa.

Figure 6 displays the pre-slip Coulomb stress zeros and the slip patch boundaries as a function of depletion pressure p for a static friction coefficient ${\mu _{st}} = 0.52$ and with all other variables as in the Example (see also Table 1). It clearly illustrates that the difference between the zeros (i.e. the potential slip patch boundaries; green curves) and the actual slip patch boundaries (blue curves) grows with increasing depletion. It also illustrates that, for a constant friction coefficient, the uncoupled approximation (indicated with a dotted blue line) performs quite well until approaching the pressure at which the two patches merge.

Figure 6. Pre-slip Coulomb stress zeros and slip patch boundaries as a function of depletion pressure $p$ for Example 1. The green curves represent the zeros ${y_i},\,i = 1, \ldots ,4,$ and the blue curves the boundaries ${\ {^\tilde{y}}_i}$ for a constant static friction coefficient ${\mu _{st}} = 0.52$ . The blue dotted curves represent the uncoupled approximation. Intersections of these curves with the vertical dashed line at a depletion pressure of -25 MPa correspond to the pre-slip Coulomb stress zeros (green horizontal lines) and slip patch boundaries (blue horizontal lines) in Figures 2 and 5. The red curves represent the boundaries ${\ {^\tilde{y}}_i}$ in case of slip weakening with static and dynamic friction coefficients ${\mu _{st}} = 0.52$ , ${\mu _{dyn}} = 0.20$ and a critical slip distance ${\delta _c} = 0.02$ m. The vertical dotted red line indicates the nucleation pressure ${p^*} = - 17.44$ MPa at which seismic slip occurs. The orange lines also represent boundaries ${\ {^\tilde{y}}_i}$ in case of slip weakening but now with a dynamic friction coefficient ${\mu _{dyn}} = 0.40$ while keeping all other parameter values the same. The corresponding nucleation pressure is indicated with a dotted orange line at ${p^*} = - 21.38$ MPa. For both of the slip-weakening cases, instability resulting in seismicity occurs at the bottom of the top patch, i.e. at boundary ${\ {^\tilde{y}}_3}$ , as shown in detail in the circular inset for the case where ${\mu _{dyn}} = 0.20$ .

Slip-weakening friction

Loss of stability

An early model to explain the occurrence of seismic events is one in which the friction coefficient drops from its static value ${\mu _{st}}$ to a lower dynamic value ${\mu _{dyn}}$ as a linear function of the ratio ${{|\delta |}}\over{{{\delta _c}}}$ , where ${\delta _c}$ is the critical slip distance. The corresponding pre-slip Coulomb stress is defined as

(40) $$\matrix{ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{{{\textstyle \sum} _C}(y,\delta (y,p),p)} { = {{\textstyle \sum} _\parallel }(y,p) - {{\textstyle \sum} _{sl}}(y,\delta (y,p),p)} \cr \, &{ { = {{\textstyle \sum} _\parallel }(y,p) - {{\textstyle \sum} ^\prime_ \bot }(y,p) \cdot \left\{ {\matrix{ {{\mu _{st}} - ({\mu _{st}} - {\mu _{dyn}}){{|\delta (y,p)|} \over {{\delta _c}}}}&{{\rm{if}}}&{|\delta (y,p)| \le {\delta _c}} \cr \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \quad\quad\quad\quad\quad\quad\quad{{\mu _{dyn}}} &{{\rm{if}}} &{|\delta (y,p)| \gt \, {\delta _c}} \cr } } \right.,}}\cr }$$

where ${\mu _{st}}$ , ${\mu _{dyn}}$ and ${\delta _c}$ may all be functions of y, and where we indicated an explicit dependence of the various stresses and the slip on the incremental pressure $p(t)$ rather than on time t directly because we consider a quasi-static situation without any explicitly time-dependent parameters.

To compute the slip patch boundaries, we can now use the same semi-analytical and numerical integration approaches as for the constant friction case, although with an iterative treatment of the dependency of ${{\textstyle \sum} _C}$ on $\delta $ . The red curve in Fig. 6 depicts the slip patch boundaries ${\tilde \!\!\!\!\,y_i},i = 1, \ldots ,4$ , as a function of p for the Example with ${\mu _{st}} = 0.52$ , ${\mu _{dyn}} = 0.20$ and ${\delta _c} = 0.02$ m. It can be seen that the patches grow in size much more rapidly than for the case with the same static friction coefficient without slip weakening (blue curve). A decrease in pressure to $p = - 17.44$ MPa leads to an unbounded value of the derivative ${\partial {\tilde{y}_3} \over {\partial p}}$ at the lower boundary of the top patch, as shown in detail in the circular inset. Quasi-static aseismic slip is now impossible for further depletion which implies that a seismic event will occur. A similar loss of stability is depicted by the orange curves which correspond to slip-weakening friction with a higher dynamic friction coefficient ( ${\mu _{dyn}} = 0.40$ ) while keeping all other parameter values the same. The less aggressive drop in friction value results in a delayed onset of seismicity at $p = - 21.38$ MPa.

A detailed analysis of the onset of seismicity under slip-weakening conditions was made by Uenishi & Rice (2003), and elaborate finite element studies using a slip-weakening model for faults in the Groningen gas field were reported by Van den Bogert (2018) and Buijze et al. (2019). In the latter two studies, the onset of seismic slip is followed by a dynamic rupture analysis which shows that the seismic slip may result in merging of the two slip patches. As demonstrated by Van den Bogert (2018) and Buijze et al. (2019), the slip in the merged patches usually stays inside the reservoir, but dynamic effects may sometimes cause the slip to propagate beyond the reservoir boundaries, in particular into the underburden.

Before the start of depletion, the pre-slip Coulomb stress just below the bottom of the reservoir (i.e. just below $y = - 150$ m) is governed by the static friction coefficient ${\mu _{st}}$ and has magnitude $\sigma _C^0 = \sigma _\parallel ^0 - \sigma _{sl}^0 = 8.57 - 15.48 = - 6.91$ MPa, where we used $\sigma _{sl}^0 = \kappa - \mu (\sigma _ \bot ^0 + \beta {p^0})$ with $\kappa = 0$ , $\mu = {\mu _{st}} = 0.52$ and $\beta = 0.90$ , and where the negative Coulomb stress implies that the fault does not slip. The corresponding pre-slip shear capacity utilisation (SCU), defined as ${\rm{SCU}} = \sigma _\parallel ^0/\sigma _{sl}^0$ , is equal to 0.55, which is way below one, the threshold above which slip will occur. The post-slip SCU just below the reservoir is defined as ${\rm{SCU}} = (\sigma _\parallel ^0 + {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel })/\sigma _{sl}^0$ where $\sigma _{sl}^0$ is now governed by the slip-dependent friction coefficient ${\mu _{dyn}} \le \mu (\delta ) \le {\mu _{st}}$ . For a dynamic friction coefficient ${\mu _{dyn}} = 0.40$ , the post-slip SCU has a lower bound $\sigma _\parallel ^0/\sigma _{sl}^0 = 0.72$ if it is assumed that the slip becomes so large that the dynamic friction coefficient is indeed reached, which is likely to happen once seismic slip occurs. For ${\mu _{dyn}} = 0.20$ , we obtain a post-slip SCU with a lower bound equal to 1.44. This latter SCU, with a value way above one, implies that fault slip governed by such a low dynamic friction coefficient would continue to propagate in a run-away fashion once it would reach the underburden. Only a heterogeneity or a change in stress state could then arrest the growth of the slip patch. Note that if such a run-away situation would occur inside the reservoir boundaries, the very large negative Coulomb stress peaks at $y = \pm 150$ m (see Fig. 2) might still form a barrier for propagation to outside the reservoir if dynamic effects would be sufficiently small. Such complexities of dynamic fault slip are outside the scope of our paper and we refer to Van den Bogert (2018), Buijze et al. (2017), Buijze et al. (2019) and Buijze (2020) for examples and further analyses.

Eigen problem

An alternative to simulation-until-seismicity is to assess the stability of the a-seismically slipped fault. Uenishi & Rice (2003) performed such a quasi-static stability analysis for a fault loaded by a single peak-shaped shear stress distribution, building on earlier work by, amongst others, Dascalu et al. (2000). Uenishi & Rice (2003) developed an explicit expression for the nucleation length $\Delta {y^*}$ , defined as the slip patch length at which stable aseismic slip will cease to be possible and seismic slip will occur. The corresponding pressure is then the nucleation pressure ${p^*}$ .

In the paper by Uenishi & Rice (2003), the pre-slip Coulomb stress has a particular form while slip weakening is defined in a slightly different way than in Equation 40, that is, in terms of the slip resistance ${{\textstyle \sum} _{sl}} = {{\textstyle \sum} ^\prime_ \bot }\mu $ rather than in terms of $\mu $ . The formulation in terms of ${{\textstyle \sum} ^\prime_ \bot }\mu $ , instead of $\mu $ , makes no difference if the normal stress ${{\textstyle \sum} ^\prime_ \bot }$ is constant along the fault. However, in our case of an inclined displaced fault, the normal stresses vary along the fault, see Fig. 2 (left). More importantly, in our case the pre-slip Coulomb stress distribution is more complex, with a combined dependence on both space and pressure. To cope with these differences, we extended the stability analysis of Uenishi & Rice (2003) as described in detail in Appendices C and D. Assuming that coupling effects may be neglected, the resulting generalised eigen problem becomes

(41) $$- W(y,p)\;\dot \delta (y,p) = A\int_{{y_ - }(p)}^{{y_ + }(p)} {\frac{{\nabla \;\dot \delta (\xi ,p)}}{{\xi - y}}d\xi } {\mkern 1mu} ,$$

where dotted variables indicate differentiation with respect to p, and $W(y,p)$ is a measure of slip weakening as defined in detail in Equation C-9 in Appendix C.

Equation 41 resembles the regular eigen equation derived by Uenishi & Rice (2003) except for the presence of a variable coefficient $W(y,p)$ , instead of just $W(p)$ , which introduces an additional dependency on y. Appendix D describes a semi-analytical approximation to this generalised eigen problem in terms of Chebyshev polynomials. The eigenvalues are (weakly) pressure-dependent, and the largest eigenvalue is equal to $\Delta {y^*}$ if the slip-induced stresses ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \sigma } _\parallel }$ obey conditions 32. Unlike in the case considered by Uenishi & Rice (2003), it seems out of reach to determine the nucleation pressure ${p^*}$ by solving the eigenproblem 41 in a stand-alone fashion.

As an alternative to the full eigenvalue computation, one may approximate Equation 41 through replacing $W(y,p)$ by a spatial average $\overline W(p) = ({y_ + } - {y_ - }{)^{ - 1}}\int_{{y_ - }}^{{y_ + }} W(y,p){\kern 1pt} dy$ , an approach that was taken by Van den Bogert (2018) and Buijze et al. (2019) in conjunction with finite element simulations. After replacing $W(y,p)$ by $\overline W(p)$ , Equation 41 then becomes the regular linear eigen equation that was solved semi-analytically by Uenishi & Rice (2003). We can therefore directly use their result which can be expressed as:

(42) $$\Delta y_{U\& R}^* = 1.158{G \over {\overline W(p)(1 - \nu )}}.$$

For simulation with a small final pressure step and a tight convergence setting, we found slightly lower results for $\Delta y_{U \& R}^*$ in comparison to $\Delta y_{eig}^*$ which illustrates that the effects of the spatial dependence of W and of the more complex ‘loading stress’ are limited. Note, however, that this approximate approach only leads to a moderate computational advantage: it skips the eigenvalue computation but still requires an iterative computation of the slip distribution to determine the spatial average $\overline W$ .

It should be noted that the nucleation length based on eigenvalues, whether $\Delta y_{eig}^*$ or $\Delta y_{U \& R}^*$ , is obtained without taking into account the coupling effect between the patches. The distance between the slip patches is influenced by the relative fault throw ${t_{f}}/h$ , that is, the fault throw normalised with respect to reservoir height. For ${t_{f}}/h$ getting closer to one, the distance becomes smaller and thus the coupling effect larger. An example of this effect will be discussed below.

Termination criteria

Simulation up to nucleation can be performed through solving the inverse Equation 30 in combination with conditions 32 while iteratively taking care of the slip dependency of the Coulomb stress, for a gradually decreasing pressure p. Close to nucleation it requires very small pressure steps and tight convergence tolerances to avoid overshooting the critical slip length $\Delta {y^*}$ and the associated nucleation pressure ${p^*}$ . Termination of the simulation can then be governed through numerically computing the derivative ${{\partial p}} \over {{\partial \Delta y}}$ for each of the slip patches and checking if it drops below a predefined tolerance close to zero; see Fig. C-1 (right) in Appendix C.

An alternative approach is to simulate the slip development and, when approaching nucleation, concurrently solve the eigen system at each pressure step. Whenever the simulated critical slip length $\Delta y_{sim}^*$ in one of the patches becomes equal to the nucleation length $\Delta y_{eig}^*$ for that patch, the simulation is terminated; see Fig. C-1 (middle) in Appendix C. However, the eigenproblem does not take into account coupling effects and may therefore overpredict $\Delta {y^*}$ and $|{p^*}|$ , a situation that will typically occur for large relative fault offsets ${t_{f}}/h$ . In practice, termination of the simulation can therefore be chosen to occur when either the numerical derivative drops below a value close to zero or when $\Delta y_{sim}^*$ overshoots $\Delta y_{eig}^*$ . Alternatively, one can simply simulate until convergence failure and accept the last converged results, but this easily leads to an over-prediction of the nucleation length unless very small pressure steps and very tight convergence control are used.

In case the slip-induced reduction in fault friction becomes so large that ${\mu _{dyn}}$ is locally reached, a further reduction will no longer take place in that location. If this happens in one patch before the critical slip length $\Delta {y^*}$ is reached, that is, before nucleation takes place, continued patch growth will occur which may now reach values beyond $\Delta {y^*}$ without leading to nucleation in that patch, although seismicity may still be triggered by nucleation in the other patch. If $\mu $ drops to ${\mu _{dyn}}$ in both patches before $\Delta {y^*}$ is reached, nucleation will be avoided altogether (Van den Bogert, 2018; Buijze, 2020). With continuing depletion, non-seismic merging of the patches may then occur. We refer to the report by Van den Bogert (2018) for a detailed analysis and a taxonomy of the various possible combinations of merging and slip.

Seismic moment

For a given nucleation pressure ${p^*}$ , Equation 39 provides an upper bound for the seismic moment per unit strike length under the assumption that slip remains confined to the reservoir. A tighter upper bound, under the same assumption, can be obtained by starting from the slip configuration ${\delta ^*}(y)$ as follows. In line with the quasi-steady-state approach that was used to simulate depletion, we can assume that the post-event configuration is in a quasi-steady state again and corresponds to a post-seismic slip distribution ${\delta _{ps}}(y)$ with post-seismic slip patch boundaries ${\widetilde y_{ \pm ,ps}}$ in a fault with a post-seismic (residual) static friction coefficient ${\mu _{ps}}$ that is equal to the former dynamic friction coefficient ${\mu _{dyn}}$ . Thus, we assume that after the onset of seismicity when $\Delta y = \Delta {y^*}$ , additional slip occurs under dynamic conditions leading to a reduction of the friction coefficient to its residual value. The approximate seismic moment can then be expressed as

(43) $${{M_s} = G\left( {\int_{{\ {{{\tilde{y}}}}_{1,ps}}}^{{\ {{{\tilde{y}}}}_{2,ps}}} {\delta _{ps}}(y)\,dy + \int_{{\ {{{\tilde{y}}}}_{3,ps}}}^{{\ {{{\tilde{y}}}}_{4,ps}}} {\delta _{ps}}(y)\,dy - \int_{\ {{{\tilde{y}}}\,}_{\!1}^*}^{\,{{\tilde{y}}}_2^*} {\delta ^*}(y)\,dy - \int_{\,{{\tilde{y}}}_3^*}^{\,\,{{\tilde{y}}}_4^*} {\delta ^*}(y)\,dy} \right)},$$

where the last two terms constitute the moment released during aseismic slip. The post-seismic slip configuration and the corresponding boundaries can be approximated by simulating aseismic slip with a constant friction coefficient ${\mu _{ps}}$ up to ${p^*}$ .

In a 3D fault configuration, the slip patches will have a finite size in the strike direction. The shape of the patches in that direction is unknown and requires a 3D analysis to be determined, especially in view of the coupling effect between the patches, the possibility of patch merging before or after nucleation, and the limits on vertical propagation beyond the reservoir boundaries in case of depletion; see also Weng & Ampuero (2019). The results from Equation 43 are therefore only relevant to illustrate the dependence of ${M_s}$ in a 2D setting. Moreover, in reality, dynamic effects may strongly influence the post-nucleation stress distribution while also more realistic friction behaviour will probably strongly influence the results. An example of ${M_s}$ versus relative fault throw ${t_{f}}/h$ will be discussed in the next section.

Computational aspects and example (continued)

We implemented the theory for simulation of fault slip using the software package Matlab with details as described in Appendix E. To compute the results in Fig. 6, representing the slip patch boundaries as a function of gradually increasing depletion, we used a variable-step-size algorithm for the depletion pressure steps. The algorithm employs a target value for the change of $\Delta y$ in the slip patches such that strong changes in function values near merging or nucleation limits can be accurately traced.

For the Example considered in this paper and for a large number of others, practical convergence was usually achieved for 100–300 Chebyshev terms with a regularisation parameter $\eta = 0.10$ m. Typical wall clock times for simulation until the onset of nucleation or merging were in the order of seconds up to tens of seconds for the semi-analytical approach, and up to minutes for the staggered grid approach, strongly depending on the setting of the various numerical parameters for spatial discretisation, regularisation, pressure stepping and convergence tolerance. Overall, the semi-analytical integration with Chebyshev polynomials proved to be computationally more robust and efficient than the staggered grid approach.

Figure 7 displays several properties as a function of fault throw ${t_{f}}$ in the Example with ${\mu _{dyn}} = 0.30$ , while keeping all other variables fixed. The fault throw varies in 45 steps between 0 and 222.5 m which corresponds to a range of dimensionless fault throw values $0 \le {t_{f}}/h \le 0.99$ . (There are no theoretical restrictions to extending this range to values above 1; the upper limit of 0.99 is a consequence of our current implementation.)

Figure 7. Results corresponding to the Example with ${\mu _{dyn}} = 0.30$ and fault throw ${t_{f}}$ varying between 0 and 222.5 m in 45 increments (relative fault throw $0 \le {t_{f}}/h \le 0.99$ ). a) nucleation pressure ${p^*}$ ; no nucleation occurs during the first four steps. b) nucleation length $\Delta {y^*}$ of top patch; blue: $\Delta y_{eig}^*$ , green: $\Delta y_{U \& R}^*$ , red: $\Delta y_{sim}^*$ ; c) seismic moment per unit strike length; blue: upper-bound ${\hat M_s}$ according to Equation 39; red: more accurate estimate ${M_s}$ according to Equation 43; d) elapsed computing time ${t_{comp}}$ . Each bar represents the time required for simulation to nucleation plus the time for an additional simulation to compute the post-seismic slip distribution.

Sub-figure (a) displays the nucleation pressure ${p^*}$ . For the first four increments, no nucleation occurs before reaching full depletion of $ - 30$ MPa. Thereafter, nucleation occurs increasingly faster, in other words, the (negative) pressure ${p^*}$ gradually drops while the nucleation length $\Delta {y^*}$ increases.

Sub-figure (b) displays the nucleation length and illustrates that the approximation $\Delta y_{U \& R}^*$ (green triangles) is pretty close to the eigenvalue-based result $\Delta y_{eig}^*$ (blue triangles). Moreover, the latter is equal to the simulated nucleation length $\Delta y_{sim}^*$ (red triangles) for increasing values of ${t_{f}}/h$ until, rather suddenly, at ${t_{f}} = 180$ m, that is, at ${t_{f}}/h = 0.80$ , the effect of coupling starts to strongly influence the growth of the two slip patches. As a result, the eigenvalue-based result $\Delta y_{eig}^*$ is no longer a good measure of the true nucleation length which rapidly drops to smaller values. We note that for all increments in Fig. 7, nucleation occurs before merging of the slip patches; in other words, they all correspond to Mechanism 2 in the taxonomy of Van den Bogert (2018), while, also, in none of these cases the value of the friction coefficient drops to the lower limit ${\mu _{dyn}}$ before nucleation occurs.

Sub-figure (c) displays the seismic moment per unit strike length ${M_s}$ (red dots) according to Equation 43 and the upper bound ${\hat M_s}$ given by Equation 39 with $p = {p^*}$ (blue dots). They both initially increase with increasing fault throw but reach a maximum whereafter they drop to lower values again. Note, however, the caveat made in Subsection 5.3.4 about the limited relevance of the magnitude of these seismic moment values.

Sub-figure (d) displays the computing time ${t_{comp}}$ for each increment. We used 150 Chebyshev terms, a minimum allowed pressure increment $\Delta p = - 100$ Pa for the slip-weakening simulations and $\Delta p = - {10^4}$ Pa for the residual friction computations, and a relative convergence tolerance for the iterative computation of the slip boundaries of ${10^{ - 8}}$ . We ran the simulations on a lap top (Intel Core i7-6600U CPU @ 260 GHz). The corresponding computing time per increment of up to 45 s is for two simulations until the nucleation pressure: one with slip-weakening friction and another one with residual friction.

Discussion and conclusions

Our paper started off with an analytical description of depletion-induced fault stresses and then moved to a semi-analytical description of fault slip, in terms of Chebyshev polynomials, and subsequently considered the effects of three relatively simple representations of friction laws while introducing successively more numerical elements. We limited our scope to quasi-static depletion-induced aseismic slip and the onset of seismicity. We did not address slow slip at constant reservoir pressures (fault creep), dynamic effects (rupture) or intermittent seismicity (stick slip). Capturing these effects would require more complex physics and in particular more advanced friction formulations such as rate and state or mechanistic models. Also, slip is known to be governed by local fluctuations in friction properties, both at small scales (Luo & Ampuero, 2018; Lebihain et al., 2021) and at larger scales because friction properties are lithology-dependent (Hunfeld et al., 2017). More advanced friction formulations and heterogeneities in friction properties can be implemented relatively straightforward in the current formulation. However, further complexities such as heterogeneous mechanical properties and 3D faults can most likely be modelled more efficiently by adopting a next level of numerical computation in the form of the boundary element method (Crouch & Starfield, 1983).

The presence of infinitely sharp stress peaks at the internal and external reservoir-fault corners in our model will not occur in reality because of a less abrupt geometry, non-elastic material behaviour and pressure diffusion between reservoir and overburden/underburden. However, relatively sharp peaks in the fault shear stresses will persist under such more realistic circumstances which means that even for small incremental pressures some amount of (aseismic) fault slip is likely to occur. Searching for an injection or depletion threshold beyond which no slip at all will occur is therefore not very meaningful, while, moreover, numerical results will likely be mesh dependent. A more relevant feature to search for seems to be the transition from near-steady-state aseismic slip to seismicity, as represented by the nucleation pressure ${p^*}$ and the associated nucleation length $\Delta {y^*}$ . However, it is well-known that the nucleation pressure is strongly dependent on friction properties and fault geometry, while for more complex friction models, such as the rate and state model, also complex patterns of intermittent seismicity may develop; see, for example, Ampuero & Rubin (2008), Cueto-Felgueroso et al. (2017). Given these uncertainties, it seems therefore more appropriate to perform sensitivity analyses or probabilistic studies aimed at defining probability distributions for the incremental pressure that may lead to seismicity, rather than to search for a single deterministic nucleation pressure.

The first aim of our paper was to gain further insight in the poro-mechanical aspects of induced seismicity in displaced faults. We confirmed many of the findings based on numerical simulation of earlier studies, in particular those by Van den Bogert (2015, 2018), Buijze et al. (2017, 2019), Buijze (2020) and Van Wees et al. (2017). A new aspect concerned the coupling between the two slip patches that develop during induced slip in a displaced fault, which affects both the forward simulation and the eigenvalue computation. Also, we relaxed an assumption in an earlier approximate expression for slip in a frictionless fault, first obtained by Bourne & Oates (2017), and showed that the original expression underestimates the present result with a factor $2(1 - \nu )$ . The same holds for the corresponding upper bound for the induced seismic moment per unit strike length. Moreover, we showed that in case of slip-weakening friction, the maximum seismic moment per unit strike length occurs for intermediate values of the dimensionless fault throw ${t_{f}}/h$ . We derived an extension of the eigenvalue computation of Uenishi & Rice (2003) for the nucleation length which can cope with pressure-dependent and location-dependent friction. We demonstrated that the result of Uenishi & Rice (2003) is a good approximation, but that both formulations fail to predict the true nucleation length at high values of the scaled fault throw because of coupling between the slip patches.

The second aim of the paper was to clarify some of theory underlying (semi-)analytical techniques for fault modelling. To this purpose, we described the use of the nucleus-of-strain concept to obtain the stress field around an infinite row of edge dislocations representing a fault (reported in the Supporting Information). We discussed aspects of Cauchy-type singular integral equations and reviewed the analytical inversion theorem of Muskhelishvili (1953). We also addressed the logarithmic singularities and jump discontinuities in the stress field around a displaced fault, the resulting lack of Hölder continuity, and regularised expressions to circumvent these issues. We discussed the use of Chebyshev polynomials to solve the governing singular integral equations for fault slip, including coupling effects between slip patches. We provided two formulations for simulation using Chebyshev polynomials. Both use (semi-)analytical integration, but while one is also based on analytical inversion (described in Appendix C), the other one uses a numerical inversion approach (described in the Supporting Information). Moreover we implemented a formulation with numerical integration, using a staggered grid approach, and analytical inversion (briefly described in Appendix E). We conclude that semi-analytical integration with analytical inversion can provide accurate results and, in our implementation, is more efficient and more robust than the other two formulations. We modified the Uenishi & Rice (2003) approach of the eigenvalue problem that governs fault stability under slip-weakening conditions to cope with pressure-dependent and location-dependent friction, and for this, we used a direct formulation in terms of Chebyshev polynomials instead of an eigenvalue expansion.

The third aim of the paper was to investigate the scope for fast simulation of induced seismicity in multiple realisations of faulted reservoirs (Monte Carlo simulation). This requires a combination of computational speed and robustness, where the simulation of millions of fault configurations would typically require run times in the sub-second domain. Our current implementation (with typical run times in the order of seconds to tens of seconds) does not fully meet the speed requirement. More importantly, the various possible combinations of reservoir and fault geometries correspond to a variation in slip patch merging and nucleation patterns that require widely varying numerical settings to be accurately simulated without fatal convergence failures. An increase in speed may be obtained through implementing Chebyshev integration with the aid of a Fast Fourier Transform (Mason & Handscomb, 2003). Also, the approximate analysis of nucleation with an energy approach, as described by Rice & Uenishi (2010), may offer an opportunity to increase the simulation speed, while the use of numerical methods based on Chebyshev polynomials, rather than the polynomials themselves, may offer another route to speed up (Golberg, 1990; Viesca & Garagash, 2018). An increased robustness will require adaptive and robust simulation control algorithms. We conclude that the current semi-analytical formulation is not yet well suited for Monte Carlo simulation, but that with some further speed up and improved simulation control that option seems to be within reach.

The last aim of the paper was to develop a semi-analytical framework for embedded fault modelling in large-scale numerical simulation tools. We conclude that this appears to be a possibility, with the potential benefit that Muskhelishvili’s analytical inversion of Cauchy integrals may help to avoid numerical convergence problems that are well-known to sometimes frustrate the quasi-static simulation of fault slip. Moreover, an implementation with Chebyshev polynomials offers scope for adaptive error control through varying the number of terms in the expansion. Whether this offers sufficient advantages in comparison to straightforward numerical integration or boundary element techniques requires further research. We aim to test implementation in our in-house finite-volume-based poro-mechanical simulation code with embedded fracture modelling capabilities (Shokrollahzadeh Behbahani et al., 2022).