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, Reference Scholz2019; Fjaer et al., Reference Fjaer, Holt, Horsrud, Raaen and Risnes2021).

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. (Reference Jansen, Singhal and Vossepoel2019) 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. (Reference Jansen, Singhal and Vossepoel2019), similar analytical expressions were derived by Lehner (Reference Lehner2019). 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. (Reference Wu, Vilarrasa, De Simone, Saaltink and Parisio2021) and Van den Hoek & Poessé (Reference Van den Hoek and 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 (Reference Lehner2019).

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. (Reference Jansen, Singhal and Vossepoel2019) 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, Reference Bourne and Oates2017; Fjaer et al., Reference Fjaer, Holt, Horsrud, Raaen and Risnes2021).

Another, semi-analytical, method to model the elastic stresses in displaced faults was presented by Van Wees et al. (Reference Van Wees, Pluymaekers, Osinga, Fokker, Van Thienen-Visser, Orlic, Wassing, Hegen and Candela2019) who based their expressions on earlier work by Okada (Reference Okada1992). The underlying solid-mechanical formulation is closely related to those used by Jansen et al. (Reference Jansen, Singhal and Vossepoel2019) and Lehner (Reference Lehner2019) 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. (Reference Jansen, Singhal and Vossepoel2019), Lehner (Reference Lehner2019), Wu et al. (Reference Wu, Vilarrasa, De Simone, Saaltink and Parisio2021) and Van den Hoek & Poessé (Reference Van den Hoek and Poessé2021). Comparison of the stresses in 2D examples using the methods of Van Wees et al. (Reference Van Wees, Pluymaekers, Osinga, Fokker, Van Thienen-Visser, Orlic, Wassing, Hegen and Candela2019), Lehner (Reference Lehner2019) and Jansen et al. (Reference Jansen, Singhal and Vossepoel2019) 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 (Reference Hettema2020) 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 (Reference Van den Bogert2015, Reference Van den Bogert2018) and subsequently by others; see Buijze et al. (Reference Buijze, Van den Bogert, Wassing, Orlic and Ten Veen2017, Reference Buijze, Van den Bogert, Wassing and Orlic2019), Van Wees et al. (Reference Van Wees, Fokker, Van Thienen-Visser, Wassing, Osinga, Orlic, Ghouri, Buijze and Pluymaekers2017, Reference Van Wees, Osinga, Van Thienen-Visser and Fokker2018, Reference Van Wees, Pluymaekers, Osinga, Fokker, Van Thienen-Visser, Orlic, Wassing, Hegen and Candela2019) for numerical studies, and Jansen et al. (Reference Jansen, Singhal and Vossepoel2019) and Lehner (Reference Lehner2019) 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. (Reference Jansen, Singhal and Vossepoel2019).

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.

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, Reference Burgers1939; Nabarro, Reference Nabarro1952; Weertman, Reference Weertman1996; Cai & Nix, Reference Cai and Nix2016). Since the 1950s, dislocations have been applied to represent earthquake sources by, for example, Vvedenskaya (Reference Vvedenskaya1956), Keylis-Borok (Reference Keylis-Borok1956) and Steketee (Reference Steketee1958a, b). We refer to Eshelby (Reference Eshelby1973), Segall (Reference Segall2010), Udías et al. (Reference Udías, Madariaga and Buforn2014) 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, Reference Nabarro1952; Barber, Reference Barber2010; Cai & Nix, Reference Cai and Nix2016). 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 (Reference Steketee1958b), Eshelby (Reference Eshelby1973), Segall (Reference Segall2010) and Udías et al. (Reference Udías, Madariaga and Buforn2014). 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, Reference Lord Kelvin1848) and the nucleus-of-strain concept (Love, Reference Love1927) 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 (Reference Bilby and Eshelby1968). 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, Reference Cates2019). 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 (Reference Muskhelishvili1953). Somewhat more accessible treatments are given by, for example, Tricomi (Reference Tricomi1957), Weertman (Reference Weertman1996) and Estrada & Kanwal (Reference Estrada and Kanwal2000). 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 (Reference Bilby and Eshelby1968) and Rice (Reference Rice1968), with a more recent treatment being given in the book by Segall (Reference Segall2010). 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.

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 (Reference Muskhelishvili1953) 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 (Reference Muskhelishvili1953) 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 (Reference Muskhelishvili1953) 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, Reference Bilby and Eshelby1968):

(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 (Reference Muskhelishvili1953), Weertman (Reference Weertman1996), Segall (Reference Segall2010) 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, Reference Tricomi1957; Weertman, Reference Weertman1996).

Equations 30 and 32, with ${s_ - } = - {s_ + }$ , were used in Bilby & Eshelby (Reference Bilby and Eshelby1968) and Rice (Reference Rice1968) to compute the stresses in various fracture configurations modelled as arrays of dislocations, and in Rice (Reference Rice1980) and Segall (Reference Segall2010) in a similar fashion to compute stresses in and around faults caused by imposed shear stresses. In Uenishi & Rice (Reference Uenishi and Rice2003), 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, Reference Weertman1996). 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 (Reference Golberg1990), Keller & Wróbel (Reference Keller and Wróbel2016) and Viesca & Garagash (Reference Viesca and Garagash2018) 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, Reference Mason and Handscomb2003; Uenishi & Rice, Reference Uenishi and Rice2003; Segall, Reference Segall2010; Barber, Reference Barber2018). 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 (Reference Uenishi and Rice2003). Further details are provided in Appendix E, which also contains a brief comparison of the computational performance of the two approaches.

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 (Reference Bourne and Oates2017) 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. (Reference Jansen, Singhal and Vossepoel2019) 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 (Reference Bourne and Oates2017), 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 (Reference Bourne and Oates2017). 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 (Reference Bourne and Oates2017) 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., Reference Novikov, Voskov, Khait, Hajibeygi and Jansen2022). 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 (Reference Bourne and Oates2017). 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 (Reference Bourne and Oates2017) 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 (Reference Uenishi and Rice2003) 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, Reference Uenishi and Rice2003). 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 (Reference Uenishi and Rice2003), and elaborate finite element studies using a slip-weakening model for faults in the Groningen gas field were reported by Van den Bogert (Reference Van den Bogert2018) and Buijze et al. (Reference Buijze, Van den Bogert, Wassing and Orlic2019). 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 (Reference Van den Bogert2018) and Buijze et al. (Reference Buijze, Van den Bogert, Wassing and Orlic2019), 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 (Reference Van den Bogert2018), Buijze et al. (Reference Buijze, Van den Bogert, Wassing, Orlic and Ten Veen2017), Buijze et al. (Reference Buijze, Van den Bogert, Wassing and Orlic2019) and Buijze (Reference Buijze2020) 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 (Reference Uenishi and Rice2003) 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. (Reference Dascalu, Ionescu and Campillo2000). Uenishi & Rice (Reference Uenishi and Rice2003) 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 (Reference Uenishi and Rice2003), 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 (Reference Uenishi and Rice2003) 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 (Reference Uenishi and Rice2003) 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 (Reference Uenishi and Rice2003), 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 (Reference Van den Bogert2018) and Buijze et al. (Reference Buijze, Van den Bogert, Wassing and Orlic2019) 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 (Reference Uenishi and Rice2003). 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.

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 (Reference Weng and Ampuero2019). 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.