Complex interplay between stress perturbations and viscoelastic relaxation in a two-asperity fault model

We consider a plane fault with two asperities embedded in a shear zone, subject to a uniform strain rate owing to tectonic loading. After an earthquake, the static stress field is relaxed by viscoelastic deformation in the asthenosphere. We treat the fault as a discrete dynamical system with 3 degrees of freedom: the slip deficits of the asperities and the variation of their difference due to viscoelastic deformation. The evolution of the fault is described in terms of inter-seismic intervals and slip episodes, which may involve the slip of a single asperity or both. We consider the effect of stress transfers connected to earthquakes produced by neighbouring faults. The perturbation alters the slip deficits of both asperities and the stress redistribution on the fault associated with viscoelastic relaxation. The interplay between the stress perturbation and the viscoelastic relaxation significantly complicates the evolution of the fault and its seismic activity. We show that the presence of viscoelastic relaxation prevents any simple correlation between the change of Coulomb stresses on the asperities and the anticipation or delay of their failures. As an application, we study the effects of the 1999 Hector Mine, California, earthquake on the post-seismic evolution of the fault that generated the 1992 Landers, California, earthquake, which we model as a two-mode event associated with the consecutive failure of two asperities.

1 Introduction

Asperity models have long been acknowledged as an effective means to describe many aspects of fault dynamics (Lay et al., 1982; Scholz, 2002). In such models, it is assumed that the bulk of energy release during an earthquake is due to the failure of one or more regions on the fault characterized by a high static friction and a velocity-weakening dynamic friction. The stress build-up on the asperities is governed by the relative motion of tectonic plates. Earthquakes that have been ascribed to the slip of two asperities are the 1964 Alaska earthquake (Christensen and Beck, 1994); the 1992 Landers, California, earthquake (Kanamori et al., 1992); the 2004 Parkfield, California, earthquake (Twardzik et al., 2012); and the 2010 Maule, Chile, earthquake (Delouis et al., 2010).

In the framework of asperity models, a critical role is played by stress accumulation on the asperities, fault slip at the asperities and stress transfer between the asperities. Accordingly, fault dynamics can be fruitfully investigated via discrete dynamical systems whose essential components are the asperities (Ruff, 1992; Turcotte, 1997). Such an approach reduces the number of degrees of freedom required to describe the dynamics of the system, that is, the evolution of the fault (in terms of slip and stress distribution) during the seismic cycle; also, it allows the visualization of the state of the fault and to follow its evolution via a geometrical approach, by means of orbits in the phase space. Finally, a finite number of dynamic modes can be defined, each one describing a particular phase of the evolution of the fault (e.g. tectonic loading, seismic slip, after-slip). Asperity models are capable of reproducing the essential features of the seismic source, while sparing the more complicated characterization based on continuum mechanics.

In a number of recent works, modelling of different mechanical phenomena in a two-asperity fault system has been addressed, such as stress perturbations due to surrounding faults (Dragoni and Piombo, 2015) and the radiation of seismic waves (Dragoni and Santini, 2015). In these models, the fault is treated as a discrete dynamical system with four dynamic modes: a sticking mode, corresponding to stationary asperities, and three slipping modes, associated with the separate or simultaneous failure of the asperities.

In the framework of a discrete fault model, the impact of viscoelastic relaxation was first studied by Amendola and Dragoni (2013) and then further investigated by Dragoni and Lorenzano (2015), who considered a fault with two asperities of different strengths. The authors discussed the features of the seismic events predicted by the model and showed how the shape of the associated source functions is related to the sequence of dynamic modes involved. In turn, the observation of the moment rate provides an insight into the state of the system at the beginning of the event, that is, the particular stress distribution on the fault from which the earthquake takes place.

However, no fault can be considered isolated; in fact, any fault is subject to stress perturbations associated with earthquakes on neighbouring faults (Harris, 1998; Stein, 1999; Steacy et al., 2005). Whenever a fault slips, the stress field in the surrounding medium is altered. As a result, the occurrence time and the magnitude of next earthquakes may change with respect to the unperturbed condition, which is governed by tectonic loading.

The aim of the present paper is to discuss the combined effects of viscoelastic relaxation and stress perturbations on a two-asperity fault in the framework of a discrete fault model. In order to deal with such a problem, we base our work on the results achieved by Dragoni and Piombo (2015) and Dragoni and Lorenzano (2015). In the former work, the authors considered a two-asperity fault with purely elastic rheology and discussed the effect of stress perturbations due to earthquakes on neighbouring faults. The fault was treated as a discrete dynamical system whose state is described by two variables, the slip deficits of the asperities. In the latter work, viscoelastic relaxation on the fault was dealt with by adding a third state variable, the variation in the difference between the slip deficits of the asperities during inter-seismic intervals. In the present paper, we introduce stress perturbations as modelled by Dragoni and Piombo (2015) in the framework of the two-asperity fault considered by Dragoni and Lorenzano (2015). Accordingly, the present work represents a three-dimensional generalization of the model devised by Dragoni and Piombo (2015). Elastic wave radiation and additional constraints on the state of the fault are taken into account, as further developments with respect to previous works.

In the framework of the present model, seismic events generated by the fault are discriminated according to the number and sequence of slipping modes involved and the seismic moment released; these features are related to the particular state of the system at the beginning of the inter-seismic interval preceding the event. We discuss how stress perturbations affect the evolution of the fault in terms of changes in the state of the system and in the duration of the inter-seismic time, highlighting the complications arising from the ongoing post-seismic deformation process with respect to the purely elastic case considered by Dragoni and Piombo (2015). As an application, we consider the stress perturbation imposed by the 1999 Hector Mine, California, earthquake (Jónsson et al., 2002; Salichon et al., 2004) to the fault that caused the 1992 Landers, California, earthquake, which we model as a two-mode event due to the consecutive failure of two asperities and that was followed by remarkable viscoelastic relaxation (Kanamori et al., 1992; Freed and Lin, 2001). We propose a means to estimate the stress transfer from the knowledge of the relative positions and faulting styles of the two faults. As a further novelty with respect to the work presented by Dragoni and Lorenzano (2015), we show how the knowledge of the time interval elapsed after the 1999 earthquake can be used to constrain the admissible set of states that may have given rise to the 1992 event. We discuss the possible subsequent evolution of the Landers fault after the stress transfer from the Hector Mine earthquake, pointing out the main differences with respect to an unperturbed scenario.

2 The model

We consider a plane fault containing two asperities of equal areas A and different strengths that we name asperity 1 and asperity 2 (Fig. 1). The fault is enclosed in a homogeneous and isotropic shear zone behaving as a Poisson solid and is subject to a uniform strain rate owing to the relative motion of two tectonic plates, taking place at a constant rate V. In order to describe the viscoelastic relaxation in the asthenosphere, we assume a Maxwell viscoelastic behaviour with a characteristic relaxation time Θ. Finally, we characterize the seismic efficiency of the fault by means of an impedance γ. All quantities are expressed in non-dimensional form.

In the present model we do not consider aseismic slip on the fault. It has been treated in the framework of a discrete fault model by Dragoni and Lorenzano (2017), who considered a region slipping aseismically for a finite time interval and calculated the effect on the stress distribution and the subsequent evolution of the fault. Of course, if the amplitude of aseismic slip has the same order of magnitude as that of seismic slip, the fault evolution may be affected.

In accordance with the assumptions of asperity models, we ascribe the generation of earthquakes on the fault to the failure of the sole asperities, neglecting any contribution of the surrounding weaker region to the seismic moment. Also, we do not describe friction, slip and stress at every point of the fault, but only consider their average values on each asperity.

Figure 1 Sketch of the model of a plane fault with two asperities. The rectangular frame is the fault border. The state of the asperities is described by their slip deficits X(T) and Y(T), while the variable Z(T) represents the temporal variation of the difference between the slip deficits of the asperities due to viscoelastic deformation in the inter-seismic interval following an earthquake on the fault.

The fault is treated as a dynamical system with three state variables, functions of time T: the slip deficits X(T) and Y(T) of asperity 1 and 2, respectively, and the variable Z(T) representing the temporal variation of the difference Y−X between the slip deficits of the asperities during inter-seismic intervals of the fault, due to the stress redistribution associated with viscoelastic relaxation in the asthenosphere. The slip deficit of an asperity is defined as the slip that an asperity should undergo at a given instant in time in order to recover the relative displacement of tectonic plates that occurred up to that moment.

We assume the simplest form of rate-dependent friction and associate the asperities with constant static and dynamic frictions, the latter considered as the average value during slip. The static friction on asperity 2 is a fraction β of that on asperity 1 and dynamic frictions are a fraction ϵ of static frictions for both asperities. Letting f_s1 and f_d1 be the static and dynamic frictions on asperity 1, respectively, and f_s2 and f_d2 be the static and dynamic frictions on asperity 2, respectively, we have

$$\begin{array}{}\text{(1)}& \mathit{\beta }={\displaystyle \frac{f_{\mathrm{s}\mathrm{2}}}{f_{\mathrm{s}\mathrm{1}}}}={\displaystyle \frac{f_{\mathrm{d}\mathrm{2}}}{f_{\mathrm{d}\mathrm{1}}}},\mathit{ϵ}={\displaystyle \frac{f_{\mathrm{d}\mathrm{1}}}{f_{\mathrm{s}\mathrm{1}}}}={\displaystyle \frac{f_{\mathrm{d}\mathrm{2}}}{f_{\mathrm{s}\mathrm{2}}}}.\end{array}$$

We acknowledge that the values of friction after a seismic event might be different from the initial ones, even though it is probable that the change is remarkable only after several seismic cycles. We neglect this possible change, because we focus on other sources of irregularity in the seismic cycles. However, the model could easily incorporate a change in friction after each event: new values could be given to static and dynamic frictions after the event and the subsequent evolution could be calculated accordingly.

During a global stick mode, the tangential forces acting on the asperities in the slip direction are (in units of static friction on asperity 1)

$$\begin{array}{}\text{(2)}& F_{\mathrm{1}}=-X+\mathit{\alpha }Z,F_{\mathrm{2}}=-Y-\mathit{\alpha }Z.\end{array}$$

In these expressions, the terms −X and −Y represent the effect of tectonic loading and have the same sign for both asperities, whereas the terms ±αZ are the contributions of stress transfer between the asperities. In the framework of the present model, stress is transferred by one asperity to the other as a result of co-seismic slip; in the subsequent inter-seismic interval, the static stress field generated by asperity slip undergoes a certain amount of relaxation owing to viscoelasticity. The parameter α is a measure of the degree of coupling between the asperities: for smaller values of α, the stress transfer from one asperity to the other is less efficient. In the limit case α=0, the asperities are completely independent from one another and the slip of one of them does not affect the state of the other: the evolution of the asperities is thus governed by tectonic loading only. By comparison with a model based on continuum mechanics, the specific value of α can be estimated as (Dragoni and Tallarico, 2016)

$$\begin{array}{}\text{(3)}& \mathit{\alpha }={\displaystyle \frac{Avs}{\mathrm{2}\dot{e}}},\end{array}$$

where A is the area of the asperities, v is the velocity of the tectonic plates, s is the tangential traction (per unit moment) imposed on one asperity by the slip of the other and $\dot{e}$ is the tangential strain rate on the fault due to tectonic loading.

An effective way to characterize fault mechanics is provided by the concept of Coulomb stress (Stein, 1999). It is defined as the difference between the shear stress σ_t in the direction of fault slip and the static friction τ_s on the fault surface:

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_C={\mathit{\sigma }}_{\mathrm{t}}-{\mathit{\tau }}_{\mathrm{s}}.\end{array}$$

Accordingly, σ_C is negative during an inter-seismic interval and a seismic event occurs when σ_C vanishes. In our model, the presence of two asperities makes it necessary to assign a value of Coulomb stress to each of them. By definition, the Coulomb forces on asperity 1 and 2 are, respectively,

$$\begin{array}{}\text{(5)}& F_{\mathrm{1}}^{\mathrm{C}}=-F_{\mathrm{1}}-\mathrm{1},F_{\mathrm{2}}^{\mathrm{C}}=-F_{\mathrm{2}}-\mathit{\beta }.\end{array}$$

Using Eq. (2), they can be rewritten as

$$\begin{array}{}\text{(6)}& F_{\mathrm{1}}^{\mathrm{C}}=X-\mathit{\alpha }Z-\mathrm{1},F_{\mathrm{2}}^{\mathrm{C}}=Y+\mathit{\alpha }Z-\mathit{\beta }.\end{array}$$

To sum up, the system is described by the set of six parameters α, β, γ, ϵ, Θ and V, with α≥0, $\mathrm{0}<\mathit{\beta }<\mathrm{1}$, γ≥0, $\mathrm{0}<\mathit{ϵ}<\mathrm{1}$, Θ>0 and V>0. At any instant T in time, the state of the system may be univocally expressed by the tern $\left(X,Y,Z\right)$ or by one of the couples $(F_{\mathrm{1}},F_{\mathrm{2}}),(F_{\mathrm{1}}^{\mathrm{C}},F_{\mathrm{2}}^{\mathrm{C}})$.

When considering the fault dynamics during the seismic cycle, it is possible to identify four dynamic modes, each one described by a different system of autonomous ODEs (ordinary differential equations): a sticking mode (00), corresponding to stationary asperities, and three slipping modes, associated with the slip of asperity 1 alone (mode 10), the slip of asperity 2 alone (mode 01) and the simultaneous slip of the asperities (mode 11). A seismic event generally consists of n slipping modes and involves one or both the asperities.

The sticking region

The sticking region of the system is defined as the set of states in which both asperities are stationary. During a global stick phase (mode 00), the rates $\dot{X}$, $\dot{Y}$ and $\dot{Z}$ are negligible with respect to their values when the asperities are slipping; thus, the sticking region is a subset of the space XYZ.

The slip of asperity 1 occurs when

$$\begin{array}{}\text{(7)}& F_{\mathrm{1}}=-\mathrm{1},\end{array}$$

while the slip of asperity 2 takes place when

$$\begin{array}{}\text{(8)}& F_{\mathrm{2}}=-\mathit{\beta }.\end{array}$$

Combining these conditions with the expressions (2) of the forces, we obtain two planes in the XYZ space,

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }X-\mathit{\alpha }Z-\mathrm{1}=\mathrm{0},\text{(10)}& {\displaystyle }& {\displaystyle }Y+\mathit{\alpha }Z-\mathit{\beta }=\mathrm{0},\end{array}$$

which we name Π₁ and Π₂, respectively. Of course, the Coulomb forces $F_{\mathrm{1}}^{\mathrm{C}}$ and $F_{\mathrm{2}}^{\mathrm{C}}$ vanish on Π₁ and Π₂, respectively; furthermore, their gradients

$$\begin{array}{}\text{(11)}& \mathrm{\nabla }{F}_{\mathrm{1}}^{\mathrm{C}}=\left(\mathrm{1},\mathrm{0},-\mathit{\alpha }\right),\mathrm{\nabla }{F}_{\mathrm{2}}^{\mathrm{C}}=\left(\mathrm{0},\mathrm{1},\mathit{\alpha }\right)\end{array}$$

are orthogonal to Π₁ and Π₂, respectively.

We exclude overshooting during the slipping modes: accordingly, we assume X≥0 and Y≥0. As a consequence, the tangential forces on the asperities must always be in the same direction as the velocity of tectonic plates, i.e. F₁≤0 and F₂≤0. From Eq. (2), the limit cases F₁=0 and F₂=0 correspond to two planes in the XYZ space,

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }X-\mathit{\alpha }Z=\mathrm{0},\text{(13)}& {\displaystyle }& {\displaystyle }Y+\mathit{\alpha }Z=\mathrm{0},\end{array}$$

which we name Γ₁ and Γ₂, respectively.

To sum up, the sticking region of the system is the subset of the XYZ space enclosed by the planes $X=\mathrm{0},Y=\mathrm{0}$, Γ₁, Γ₂, Π₁ and Π₂: a convex hexahedron H. Its vertices are the origin $(\mathrm{0},\mathrm{0},\mathrm{0})$ and the points

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }A=\left(\mathrm{0},\mathrm{1},-{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}\right),B=\left(\mathit{\beta },\mathrm{0},{\displaystyle \frac{\mathit{\beta }}{\mathit{\alpha }}}\right),\\ \text{(14)}& {\displaystyle }& {\displaystyle }C=\left(\mathit{\beta }+\mathrm{1},\mathrm{0},{\displaystyle \frac{\mathit{\beta }}{\mathit{\alpha }}}\right),\text{(15)}& {\displaystyle }& {\displaystyle }D=\left(\mathrm{0},\mathit{\beta }+\mathrm{1},-{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}\right),E=\left(\mathrm{1},\mathrm{0},\mathrm{0}\right),F=\left(\mathrm{0},\mathit{\beta },\mathrm{0}\right).\end{array}$$

The sticking region is shown in Fig. 2 for a particular choice of the parameters α and β. Its volume can be expressed as a function of the parameters of the system as $\mathit{\beta }(\mathit{\beta }+\mathrm{1})/\mathrm{2}\mathit{\alpha }$. Accordingly, the subset of the state space corresponding to stationary asperities decreases with the degree of coupling between the asperities and with the asymmetry of the system (β→0). By definition, every orbit of mode 00 is enclosed within H and eventually reaches one of its faces AECD or BCDF, belonging to the planes Π₁ and Π₂, respectively, where an earthquake starts. In these cases, the system switches from mode 00 to mode 10 or mode 01, respectively. In the particular case in which the orbit of mode 00 reaches the edge CD, the system passes from mode 00 to mode 11.

3 Dynamic modes and slip in a seismic event

Let P₀∈H be the state of the system at the beginning of an inter-seismic interval. The specific location of P₀ inside the sticking region allows the prediction of the first slipping mode involved in the next seismic event on the fault. In fact, Dragoni and Lorenzano (2015) illustrated the existence of a transcendental surface Σ within H, expressed by the equation

$$\begin{array}{}\text{(16)}& V\mathrm{\Theta }\left[W\left({\mathit{\gamma }}_{\mathrm{1}}\right)-W\left({\mathit{\gamma }}_{\mathrm{2}}\right)\right]+Y-X+\mathrm{1}-\mathit{\beta }=\mathrm{0},\end{array}$$

where W is the Lambert function with arguments

$$\begin{array}{}\text{(17)}& {\mathit{\gamma }}_{\mathrm{1}}={\displaystyle \frac{\mathit{\alpha }Z}{V\mathrm{\Theta }}}{e}^{-\frac{\mathrm{1}-X}{V\mathrm{\Theta }}},{\mathit{\gamma }}_{\mathrm{2}}=-{\displaystyle \frac{\mathit{\alpha }Z}{V\mathrm{\Theta }}}{e}^{-\frac{\mathit{\beta }-Y}{V\mathrm{\Theta }}}.\end{array}$$

The surface Σ divides H in two subsets H₁ and H₂ (Fig. 3). The seismic event starts with mode 10 if P₀∈H₁ or with mode 01 if P₀∈H₂; in the particular case in which P₀∈Σ, the seismic event starts with mode 11.

Figure 2 The sticking region of the system, defined as the subset of the state space XYZ in which both asperities are at rest: a convex hexahedron H $(\mathit{\alpha }=\mathrm{1},\mathit{\beta }=\mathrm{1})$.

Mode 00 terminates at a point P₁ on the face AECD or BCDF of H. The number and sequence of slipping modes involved in the subsequent seismic event can be discriminated from the specific position of P₁. If we consider the face AECD (Fig. 4), the earthquake will be a one-mode event 10 if P₁ belongs to the trapezoid Q₁; it will be a two-mode event 10-01 if P₁ belongs to the segment s₁; it will be a three-mode event 10-11-01 or 10-11-10 if P₁ belongs to the trapezoid R₁, where the precise sequence must be evaluated numerically and depends on the particular combination of the parameters α, β, γ and ϵ. The remaining portion of the face would lead to overshooting. Analogous considerations can be made for the subsets Q₂, s₂ and R₂ on the face BCDF. In the particular case in which P₁ belongs to the edge CD, the earthquake will be a two-mode event 11-01. In fact, by definition, friction on asperity 2 is smaller than friction on asperity 1 ($\mathrm{0}<\mathit{\beta }<\mathrm{1}$); if the asperities start slipping simultaneously, asperity 1 is then bound to stop the first, while asperity 2 continues to slip. As a result, mode 11 is followed by mode 01 and the slip of the weaker asperity has a longer duration. The opposite would hold if asperity 2 were stronger than asperity 1 (β>1), so that the slip event resulting from states P₁∈CD would be a two-mode event 11-10.

Figure 3 The surface Σ that splits the sticking region H in two subsets H₁ (below) and H₂ (above) (α=1, β=1, VΘ=1). It allows the first slipping mode during an earthquake to be discriminated.

Figure 4 The faces AECD and BCDF of the sticking region H, where seismic events start, and their subsets (α=1, β=1, γ=1, ϵ=0.7). The events taking place on the face AECD (BCDF) start with mode 10 (01). The trapezoids Q_i correspond to events involving a single asperity; the segments s_i correspond to events associated with the consecutive slips of the asperities; the trapezoids R_i correspond to events involving the simultaneous slips of the asperities.

In addition, the knowledge of the position of P₁ allows the total amount of slip of the asperities and the seismic moment associated with the earthquake to be established. Let us consider an event made up of n dynamic modes and let $P_i=(X_i,Y_i,Z_i)$ be the state of the system at time T=T_i, when the ith mode starts $(i=\mathrm{1},\mathrm{2},\mathrm{\dots }n)$. The final slip amplitudes of asperity 1 and 2 are, respectively,

$$\begin{array}{}\text{(18)}& U_{\mathrm{1}}=X_{\mathrm{1}}-X_{n+\mathrm{1}},U_{\mathrm{2}}=Y_{\mathrm{1}}-Y_{n+\mathrm{1}}.\end{array}$$

Accordingly, the final seismic moment can be calculated as

$$\begin{array}{}\text{(19)}& M_{\mathrm{0}}=M_{\mathrm{1}}{\displaystyle \frac{U_{\mathrm{1}}+U_{\mathrm{2}}}{U}},\end{array}$$

where M₁ and U are the seismic moment and slip amplitude associated with a one-mode event 10 in the limit case γ=0, respectively, with

$$\begin{array}{}\text{(20)}& U=\mathrm{2}{\displaystyle \frac{\mathrm{1}-\mathit{ϵ}}{\mathrm{1}+\mathit{\alpha }}}.\end{array}$$

The possible values of U₁,U₂ and M₀ are summarized in Table 1: the effect of wave radiation is characterized by means of the quantity

$$\begin{array}{}\text{(21)}& \mathit{\kappa }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}+e^{-\frac{\mathit{\gamma }{T}_{\mathrm{s}}}{\mathrm{2}}}\right),\end{array}$$

where T_s is the duration of slip in a one-mode event (Dragoni and Santini, 2015).

As for the evolution of the variable Z(T) during the earthquake, it changes according to the equation $\ddot{Z}=\ddot{Y}-\ddot{X}$, since the relaxation process is negligible during the slip of the asperities.

4 Stress perturbations from neighbouring faults

We now consider the perturbations of the state of the fault caused by the coseismic slip on surrounding faults. Following Dragoni and Piombo (2015), we assume that (1) the perturbations occur during an inter-seismic interval; (2) the stress transfer takes place over a time interval negligible with respect to the duration of the inter-seismic interval; and (3) at the time of the perturbation, the state of the fault is sufficiently far from the failure conditions and the stress transfer is small enough that the onset of motion of either asperity is not achieved immediately.

Table 1 Final slip amplitudes U₁ and U₂ of asperity 1 and 2 and seismic moment M₀ during a seismic event made up of n slipping modes, as a function of the state P₁ where the event started. The entry e.n. is the abbreviation for evaluated numerically.

Let $\left(X,Y,Z\right)\in \mathbf{H}$ be the state of the fault at the time of the perturbation. Generally speaking, the system undergoes a transition to a new state

$$\begin{array}{}\text{(22)}& \left(X^{\prime },Y^{\prime },Z^{\prime }\right)=\left(X,Y,Z\right)+\left(\mathrm{\Delta }X,\mathrm{\Delta }Y,\mathrm{\Delta }Z\right).\end{array}$$

Since the stress transfer takes place over a time interval short with respect to the inter-seismic interval (assumption 2), viscoelastic relaxation is negligible during the perturbation and the rheology can be reasonably considered as purely elastic as the perturbation takes place. Accordingly, we set

$$\begin{array}{}\text{(23)}& \mathrm{\Delta }Z=\mathrm{\Delta }Y-\mathrm{\Delta }X.\end{array}$$

The change of state is then associated with a vector in the XYZ space:

$$\begin{array}{}\text{(24)}& \mathrm{\Delta }\mathit{R}=\left(\mathrm{\Delta }X,\mathrm{\Delta }Y,\mathrm{\Delta }Z\right).\end{array}$$

The components of ΔR generally have different magnitudes and may have different signs, as a consequence of the inhomogeneity of the stress field produced by an earthquake. They can be written in terms of the tangential forces ΔF₁ and ΔF₂ exerted by the perturbing source on asperity 1 and 2, respectively: from Eq. (2), we have

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{F}_{\mathrm{1}}=-\mathrm{\Delta }X+\mathit{\alpha }\mathrm{\Delta }Z=\mathit{\alpha }\mathrm{\Delta }Y-(\mathrm{1}+\mathit{\alpha })\mathrm{\Delta }X,\text{(26)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{F}_{\mathrm{2}}=-\mathrm{\Delta }Y-\mathit{\alpha }\mathrm{\Delta }Z=\mathit{\alpha }\mathrm{\Delta }X-(\mathrm{1}+\mathit{\alpha })\mathrm{\Delta }Y.\end{array}$$

Combining these expressions together, we get

$$\begin{array}{}\text{(27)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }X=-{\displaystyle \frac{\mathrm{1}+\mathit{\alpha }}{\mathrm{1}+\mathrm{2}\mathit{\alpha }}}\mathrm{\Delta }{F}_{\mathrm{1}}-{\displaystyle \frac{\mathit{\alpha }}{\mathrm{1}+\mathrm{2}\mathit{\alpha }}}\mathrm{\Delta }{F}_{\mathrm{2}},\text{(28)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }Y=-{\displaystyle \frac{\mathit{\alpha }}{\mathrm{1}+\mathrm{2}\mathit{\alpha }}}\mathrm{\Delta }{F}_{\mathrm{1}}-{\displaystyle \frac{\mathrm{1}+\mathit{\alpha }}{\mathrm{1}+\mathrm{2}\mathit{\alpha }}}\mathrm{\Delta }{F}_{\mathrm{2}},\text{(29)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }Z={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathrm{2}\mathit{\alpha }}}\left(\mathrm{\Delta }{F}_{\mathrm{1}}-\mathrm{\Delta }{F}_{\mathrm{2}}\right).\end{array}$$

We conclude that the variations in tangential stress alter the orbit of the system.

The components of ΔR can also be related to the orientation of the vector in the state space. With reference to Fig. 5, we have

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{\Delta }X=\mathrm{\Delta }R\cos \mathit{\delta }\cos \mathit{\theta },\mathrm{\Delta }Y=\mathrm{\Delta }R\cos \mathit{\delta }\sin \mathit{\theta },\\ \text{(30)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }Z=\mathrm{\Delta }R\sin \mathit{\delta }.\end{array}$$

Introducing the assumption (17), the angle δ may be expressed in terms of the angle θ as

$$\begin{array}{}\text{(31)}& \mathit{\delta }=\arctan \left(\sin \mathit{\theta }-\cos \mathit{\theta }\right).\end{array}$$

In writing Eq. (19), we took into account that

$$\begin{array}{}\text{(32)}& \mathit{\delta }\ne {\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}},{\displaystyle \frac{\mathrm{3}\mathit{\pi }}{\mathrm{2}}}\end{array}$$

or it would result in

$$\begin{array}{}\text{(33)}& \mathrm{\Delta }Z=\pm \mathrm{\Delta }R,\mathrm{\Delta }X=\mathrm{\Delta }Y=\mathrm{0},\end{array}$$

which is a meaningless circumstance. From Eq. (19), the tangential forces (19)–(19) can be rewritten as

$$\begin{array}{}\text{(34)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{F}_{\mathrm{1}}={\displaystyle \frac{\mathit{\alpha }\sin \mathit{\theta }-(\mathrm{1}+\mathit{\alpha })\cos \mathit{\theta }}{\sqrt{\mathrm{2}-\sin \mathrm{2}\mathit{\theta }}}}\mathrm{\Delta }R,\text{(35)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{F}_{\mathrm{2}}={\displaystyle \frac{\mathit{\alpha }\cos \mathit{\theta }-(\mathrm{1}+\mathit{\alpha })\sin \mathit{\theta }}{\sqrt{\mathrm{2}-\sin \mathrm{2}\mathit{\theta }}}}\mathrm{\Delta }R.\end{array}$$

Following the variations in normal stress, the static and dynamic frictions on each asperity are altered. Letting $f_{\mathrm{s}\mathrm{1}}^{\prime }$ and $f_{\mathrm{s}\mathrm{2}}^{\prime }$ be the new static frictions on asperity 1 and 2, respectively, we define

$$\begin{array}{}\text{(36)}& {\mathit{\beta }}_{\mathrm{1}}={\displaystyle \frac{f_{\mathrm{s}\mathrm{1}}^{\prime }}{f_{\mathrm{s}\mathrm{1}}}},{\mathit{\beta }}_{\mathrm{2}}={\displaystyle \frac{f_{\mathrm{s}\mathrm{2}}^{\prime }}{f_{\mathrm{s}\mathrm{1}}}}.\end{array}$$

The changes in static frictions are then

$$\begin{array}{}\text{(37)}& \mathrm{\Delta }{\mathit{\beta }}_{\mathrm{1}}={\mathit{\beta }}_{\mathrm{1}}-\mathrm{1},\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{2}}={\mathit{\beta }}_{\mathrm{2}}-\mathit{\beta }\end{array}$$

on asperity 1 and 2, respectively.

Figure 5 The vector ΔR and its orientation in the state space XYZ, characterizing the stress perturbation imposed on the system by earthquakes produced by neighbouring faults.

Since the stress perturbation does not alter the friction coefficients of rocks, it is reasonable to assume that the ratio ϵ between dynamic and static friction is unchanged on both asperities. Therefore, letting $f_{\mathrm{d}\mathrm{1}}^{\prime }$ and $f_{\mathrm{d}\mathrm{2}}^{\prime }$ be the new dynamic frictions on asperity 1 and 2, respectively, we have

$$\begin{array}{}\text{(38)}& {\displaystyle \frac{f_{\mathrm{d}\mathrm{1}}^{\prime }}{f_{\mathrm{s}\mathrm{1}}}}=\mathit{ϵ}{\displaystyle \frac{f_{\mathrm{s}\mathrm{1}}^{\prime }}{f_{\mathrm{s}\mathrm{1}}}}=\mathit{ϵ}{\mathit{\beta }}_{\mathrm{1}},{\displaystyle \frac{f_{\mathrm{d}\mathrm{2}}^{\prime }}{f_{\mathrm{s}\mathrm{1}}}}=\mathit{ϵ}{\displaystyle \frac{f_{\mathrm{s}\mathrm{2}}^{\prime }}{f_{\mathrm{s}\mathrm{1}}}}=\mathit{ϵ}{\mathit{\beta }}_{\mathrm{2}}.\end{array}$$

The consequent changes in dynamic frictions are ϵΔβ₁ and ϵΔβ₂ on asperity 1 and 2, respectively.

4.1 Effects of the perturbation

The stress transfer resulting from earthquakes on neighbouring faults alters several parameters of the model. A first remarkable change concerns the strength of the asperities. After the perturbation, we can define a new ratio

$$\begin{array}{}\text{(39)}& {\mathit{\beta }}^{\prime }={\displaystyle \frac{f_{\mathrm{s}\mathrm{2}}^{\prime }}{f_{\mathrm{s}\mathrm{1}}^{\prime }}}={\displaystyle \frac{f_{\mathrm{d}\mathrm{2}}^{\prime }}{f_{\mathrm{d}\mathrm{1}}^{\prime }}}={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{2}}}{{\mathit{\beta }}_{\mathrm{1}}}},\end{array}$$

which differs from the original value of β given in Eq. (1). Moreover, the stress transfer may be so intense that the weaker asperity may become the stronger one: that is, it may result in ${\mathit{\beta }}^{\prime }>\mathrm{1}$.

The variations in static frictions entail different conditions for the onset of motion of the asperities. Taking Eq. (22) into account, Eqs. (7) and (8) become, respectively,

$$\begin{array}{}\text{(40)}& F_{\mathrm{1}}=-{\mathit{\beta }}_{\mathrm{1}},F_{\mathrm{2}}=-{\mathit{\beta }}_{\mathrm{2}}.\end{array}$$

By combination with Eq. (2), these conditions define the planes

$$\begin{array}{}\text{(41)}& {\displaystyle }& {\displaystyle }X-\mathit{\alpha }Z-{\mathit{\beta }}_{\mathrm{1}}=\mathrm{0},\text{(42)}& {\displaystyle }& {\displaystyle }Y+\mathit{\alpha }Z-{\mathit{\beta }}_{\mathrm{2}}=\mathrm{0}\end{array}$$

that we call ${\mathrm{\Pi }}_{\mathrm{1}}^{\prime }$ and ${\mathrm{\Pi }}_{\mathrm{2}}^{\prime }$, respectively. Conversely, the planes Γ₁ and Γ₂ are not affected by the stress perturbation, since they do not depend on frictions. We conclude that changes in normal stress modify the sticking region of the system, describing a new hexahedron H^′ in the state space. The coordinates of its vertices are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{A}^{\prime }=\left(\mathrm{0},{\mathit{\beta }}_{\mathrm{1}},-{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{1}}}{\mathit{\alpha }}}\right),B^{\prime }=\left({\mathit{\beta }}_{\mathrm{2}},\mathrm{0},{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{2}}}{\mathit{\alpha }}}\right),\\ \text{(43)}& {\displaystyle }& {\displaystyle }{C}^{\prime }=\left({\mathit{\beta }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{2}},\mathrm{0},{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{2}}}{\mathit{\alpha }}}\right),{\displaystyle }& {\displaystyle }{D}^{\prime }=\left(\mathrm{0},{\mathit{\beta }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{2}},-{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{1}}}{\mathit{\alpha }}}\right),E^{\prime }=\left({\mathit{\beta }}_{\mathrm{1}},\mathrm{0},\mathrm{0}\right),\\ \text{(44)}& {\displaystyle }& {\displaystyle }{F}^{\prime }=\left(\mathrm{0},{\mathit{\beta }}_{\mathrm{2}},\mathrm{0}\right).\end{array}$$

The volume of H^′ is ${\mathit{\beta }}_{\mathrm{1}}{\mathit{\beta }}_{\mathrm{2}}({\mathit{\beta }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{2}})/\mathrm{2}\mathit{\alpha }$: thus, the set of states corresponding to stationary asperities is enlarged or reduced, depending on how normal stresses on the asperities are modified.

Following the changes in static frictions, the surface Σ in Eq. (10) is replaced by a new surface Σ^′ expressed by

$$\begin{array}{}\text{(45)}& V\mathrm{\Theta }\left[W\left({\mathit{\gamma }}_{\mathrm{1}}^{\prime }\right)-W\left({\mathit{\gamma }}_{\mathrm{2}}^{\prime }\right)\right]+Y-X+{\mathit{\beta }}_{\mathrm{1}}-{\mathit{\beta }}_{\mathrm{2}}=\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(46)}& {\mathit{\gamma }}_{\mathrm{1}}^{\prime }={\displaystyle \frac{\mathit{\alpha }Z}{V\mathrm{\Theta }}}{e}^{-\frac{{\mathit{\beta }}_{\mathrm{1}}-X}{V\mathrm{\Theta }}},{\mathit{\gamma }}_{\mathrm{2}}^{\prime }=-{\displaystyle \frac{\mathit{\alpha }Z}{V\mathrm{\Theta }}}{e}^{-\frac{{\mathit{\beta }}_{\mathrm{2}}-Y}{V\mathrm{\Theta }}}.\end{array}$$

As a result, the sticking region H^′ is split in two subsets ${\mathbf{H}}_{\mathrm{1}}^{\prime }$ and ${\mathbf{H}}_{\mathrm{2}}^{\prime }$; furthermore, its faces $A^{\prime }{E}^{\prime }{C}^{\prime }{D}^{\prime }$ and $B^{\prime }{C}^{\prime }{D}^{\prime }{F}^{\prime }$ are divided into subsets ${\mathbf{Q}}_{\mathrm{1}}^{\prime },{\mathbf{s}}_{\mathrm{1}}^{\prime }$, ${\mathbf{R}}_{\mathrm{1}}^{\prime }$ and ${\mathbf{Q}}_{\mathrm{2}}^{\prime },{\mathbf{s}}_{\mathrm{2}}^{\prime }$, ${\mathbf{R}}_{\mathrm{2}}^{\prime }$, respectively.

As a consequence of the changes in dynamic frictions, the amount of slip that asperities undergo during a seismic event is modified. In turn, the perturbation alters the seismic moment associated with an earthquake. The variations in the final slip amplitudes U₁ and U₂ of asperity 1 and asperity 2, respectively, and in the final seismic moment M₀ associated with the different seismic events predicted by the model are listed in Table 2.

4.1.1 Changes in Coulomb forces

The variations in tangential stresses and static frictions discussed so far entail a change in the Coulomb forces assigned to the asperities. Combining Eq. (5) with Eqs. (19) and (19), these changes are given by

$$\begin{array}{}\text{(47)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{F}_{\mathrm{1}}^{\mathrm{C}}=-\mathrm{\Delta }{F}_{\mathrm{1}}-\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{1}}=(\mathrm{1}+\mathit{\alpha })\mathrm{\Delta }X-\mathit{\alpha }\mathrm{\Delta }Y-\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{1}},\text{(48)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{F}_{\mathrm{2}}^{\mathrm{C}}=-\mathrm{\Delta }{F}_{\mathrm{2}}-\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{2}}=(\mathrm{1}+\mathit{\alpha })\mathrm{\Delta }Y-\mathit{\alpha }\mathrm{\Delta }X-\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{2}}\end{array}$$

or, exploiting Eqs. (22) and (22),

$$\begin{array}{}\text{(49)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{F}_{\mathrm{1}}^{\mathrm{C}}={\displaystyle \frac{(\mathrm{1}+\mathit{\alpha })\cos \mathit{\theta }-\mathit{\alpha }\sin \mathit{\theta }}{\sqrt{\mathrm{2}-\sin \mathrm{2}\mathit{\theta }}}}\mathrm{\Delta }R-\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{1}},\text{(50)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{F}_{\mathrm{2}}^{\mathrm{C}}={\displaystyle \frac{(\mathrm{1}+\mathit{\alpha })\sin \mathit{\theta }-\mathit{\alpha }\cos \mathit{\theta }}{\sqrt{\mathrm{2}-\sin \mathrm{2}\mathit{\theta }}}}\mathrm{\Delta }R-\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{2}}.\end{array}$$

The sign of $\mathrm{\Delta }{F}_i^{\mathrm{C}}$ $(i=\mathrm{1},\mathrm{2})$ determines whether the perturbation brings an asperity closer to or farther from the failure; specifically, positive variations entail that slip is favoured, and vice versa. Equations (29) and (29) clearly point out that this effect is regulated by the orientation of the vector ΔR in the state space. Bearing in mind the observations made in Sect. 2, we find that $\mathrm{\Delta }{F}_{\mathrm{1}}^{\mathrm{C}}$ is at its maximum when ΔR is perpendicular to plane Π₁ and points toward it; it vanishes when ΔR is parallel to plane Π₁, and it is at its minimum when ΔR is perpendicular to plane Π₁ and points away from it. Analogous considerations can be made for $\mathrm{\Delta }{F}_{\mathrm{2}}^{\mathrm{C}}$.

On the whole, the effect of the stress perturbation can be discussed in terms of the quantity

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }{F}^{\mathrm{C}}& {\displaystyle }=\mathrm{\Delta }{F}_{\mathrm{2}}^{\mathrm{C}}-\mathrm{\Delta }{F}_{\mathrm{1}}^{\mathrm{C}}\\ \text{(51)}& {\displaystyle }& {\displaystyle }=(\mathrm{1}+\mathrm{2}\mathit{\alpha })\left(\mathrm{\Delta }Y-\mathrm{\Delta }X\right)+\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{1}}-\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{2}}.\end{array}$$

Let us assume that the system is at a certain state $\left(X,Y,Z\right)\in {\mathbf{H}}_{\mathrm{1}}$ before the perturbation; accordingly, the next seismic event on the fault will start with the failure of asperity 1. If ΔF^C>0, the perturbation favours the slip of asperity 2 more than the slip of asperity 1: therefore, the system is brought to a state closer to the condition for the simultaneous failure of the asperities and thus to the Σ surface. On the contrary, perturbations for which ΔF^C<0 take the system farther from the Σ surface. The opposite holds for an unperturbed state $\left(X,Y,Z\right)\in {\mathbf{H}}_{\mathrm{2}}$.

Table 2 Changes in the final slip amplitudes U₁ and U₂ of asperity 1 and 2 and in the seismic moment M₀ associated with the different seismic events predicted by the model, after a stress perturbation from neighbouring faults. The entry e.n. is the abbreviation for evaluated numerically.

4.1.2 Changes in the duration of the inter-seismic interval

As already stated, stress perturbations can anticipate or delay the occurrence of an earthquake produced by a certain asperity. We now quantify this effect in terms of the variation in the duration of the inter-seismic interval. Generally speaking, the perturbation vector ΔR may cross the Σ surface and thus bring the system from an unperturbed state within H₁ (H₂) to a perturbed state within ${\mathbf{H}}_{\mathrm{2}}^{\prime }$ (${\mathbf{H}}_{\mathrm{1}}^{\prime }$). For the sake of simplicity, we consider here only the particular case in which the perturbation vector ΔR does not cross the Σ surface. An example of a more general case will be shown in Sect. 5 for a real fault.

Let us first focus on the case in which the unperturbed state $\left(X,Y,Z\right)\in {\mathbf{H}}_{\mathrm{1}}$. The time required by the orbit of the system to reach plane Π₁, triggering the failure of asperity 1, was calculated by Amendola and Dragoni (2013) as

$$\begin{array}{}\text{(52)}& T_{\mathrm{1}}=\mathrm{\Theta }W\left({\mathit{\gamma }}_{\mathrm{1}}\right)+{\displaystyle \frac{\mathrm{1}-X}{V}},\end{array}$$

with γ₁ given in Eq. (11). If the stress perturbation brings the system to a state $\left(X^{\prime },Y^{\prime },Z^{\prime }\right)\in {\mathbf{H}}_{\mathrm{1}}^{\prime }$ and the static friction on asperity 1 to β₁, the time required by the orbit to reach plane ${\mathrm{\Pi }}_{\mathrm{1}}^{\prime }$ is

$$\begin{array}{}\text{(53)}& T_{\mathrm{1}}^{\prime }=\mathrm{\Theta }W\left({\mathit{\gamma }}_{\mathrm{1}}^{\prime }\right)+{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{1}}-X^{\prime }}{V}},\end{array}$$

with ${\mathit{\gamma }}_{\mathrm{1}}^{\prime }$ given in Eq. (28). The difference between the two times is

$$\begin{array}{}\text{(54)}& \mathrm{\Delta }{T}_{\mathrm{1}}=T_{\mathrm{1}}^{\prime }-T_{\mathrm{1}}=\mathrm{\Theta }\left[W\left({\mathit{\gamma }}_{\mathrm{1}}^{\prime }\right)-W\left({\mathit{\gamma }}_{\mathrm{1}}\right)\right]-{\displaystyle \frac{\mathrm{\Delta }{F}_{\mathrm{1}}^{\mathrm{C}}+\mathit{\alpha }\mathrm{\Delta }Z}{V}},\end{array}$$

where Eq. (29) has been employed.

If instead $\left(X,Y,Z\right)\in {\mathbf{H}}_{\mathrm{2}}$, the time required by the orbit of the system to reach plane Π₂, triggering the failure of asperity 2, is given by (Amendola and Dragoni, 2013)

$$\begin{array}{}\text{(55)}& T_{\mathrm{2}}=\mathrm{\Theta }W\left({\mathit{\gamma }}_{\mathrm{2}}\right)+{\displaystyle \frac{\mathit{\beta }-Y}{V}},\end{array}$$

with γ₂ given in Eq. (11). If the stress perturbation takes the system to a state $\left(X^{\prime },Y^{\prime },Z^{\prime }\right)\in {\mathbf{H}}_{\mathrm{2}}^{\prime }$ and the static friction on asperity 2 to β₂, the time required to reach plane ${\mathrm{\Pi }}_{\mathrm{2}}^{\prime }$ is

$$\begin{array}{}\text{(56)}& T_{\mathrm{2}}^{\prime }=\mathrm{\Theta }W\left({\mathit{\gamma }}_{\mathrm{2}}^{\prime }\right)+{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{2}}-Y^{\prime }}{V}},\end{array}$$

with ${\mathit{\gamma }}_{\mathrm{2}}^{\prime }$ given in Eq. (28). The difference between the two times is

$$\begin{array}{}\text{(57)}& \mathrm{\Delta }{T}_{\mathrm{2}}=T_{\mathrm{2}}^{\prime }-T_{\mathrm{2}}=\mathrm{\Theta }\left[W\left({\mathit{\gamma }}_{\mathrm{2}}^{\prime }\right)-W\left({\mathit{\gamma }}_{\mathrm{2}}\right)\right]-{\displaystyle \frac{\mathrm{\Delta }{F}_{\mathrm{2}}^{\mathrm{C}}-\mathit{\alpha }\mathrm{\Delta }Z}{V}},\end{array}$$

where Eq. (29) has been employed. Positive values of ΔT₁ and ΔT₂ correspond to a delay in the occurrence of an earthquake on asperity 1 and 2, respectively, and vice versa.

4.1.3 Discussion

According to the model, rock rheology plays a critical role in the response to stress perturbations. In the case of purely elastic coupling between the asperities, Dragoni and Piombo (2015) showed that the changes in the duration of the inter-seismic interval prior to the failure of asperity 1 and 2 are, respectively,

$$\begin{array}{}\text{(58)}& \mathrm{\Delta }{T}_{\mathrm{1}}=-{\displaystyle \frac{\mathrm{\Delta }{F}_{\mathrm{1}}^{\mathrm{C}}}{V}},\mathrm{\Delta }{T}_{\mathrm{2}}=-{\displaystyle \frac{\mathrm{\Delta }{F}_{\mathrm{2}}^{\mathrm{C}}}{V}}.\end{array}$$

Accordingly, an increase in the Coulomb force associated with a given asperity ($\mathrm{\Delta }{F}_i^{\mathrm{C}}>\mathrm{0}$) directly yields the anticipation of the slip of that asperity, and vice versa. What is more, the variation in the duration of the inter-seismic interval is proportional to the change in the Coulomb force associated with the asperity.

Conversely, in the viscoelastic case there is no straightforward connection between the sign of $\mathrm{\Delta }{F}_i^{\mathrm{C}}$ and the anticipation or delay of an earthquake on the associated asperity. In fact, the expressions (31) and (34) obtained for ΔT₁ and ΔT₂ in the previous section indicate that the net effect depends in a non-trivial way on the particular state of the fault at the time of the stress perturbation and right after it. This result points out the complex interplay between the post-seismic evolution of a fault in the presence of viscoelastic relaxation and the stress transfer from neighbouring faults.

5 An application: perturbation of the 1992 Landers fault by the 1999 Hector Mine earthquake

We study the effects of the 16 October 1999 M_w 7.1 Hector Mine, California, earthquake on the post-seismic evolution of the fault that generated the 28 June 1992 M_w 7.3 Landers, California, earthquake. The geometry of the two faults is shown in Fig. 6.

The 1992 Landers earthquake was a right-lateral strike-slip event that can be approximated as the result of the slip of two coplanar asperities (Kanamori et al., 1992): a northern one (asperity 1) and a southern one (asperity 2), with average slips u₁=6 m and u₂=3 m, respectively. Following Dragoni and Tallarico (2016), we assume a common area A=300 km² for both the asperities. We place the centres of asperity 1 and asperity 2 at 34.46^∘ N, 116.52^∘ W and 34.20^∘ N, 116.44^∘ W, respectively, with a common depth of 8 km. The earthquake started with the failure of asperity 2, followed by the failure of asperity 1. We characterize the event by strike, dip and rake angles of 345, 85 and 180^∘, respectively, an average of the values provided by Kanamori et al. (1992) for the two phases of the earthquake.

The 1992 event was followed by remarkable post-seismic deformation, which can be interpreted as the result of several processes. For the sake of the present application, we assume viscoelastic relaxation as the most significant mechanism. We assign a viscosity $\mathit{\eta }=\mathrm{5}\times {\mathrm{10}}^{\mathrm{18}}\mathrm{Pa}\mathrm{s}$ to the lower crust at Landers, averaging the estimates provided by Deng et al. (1998), Pollitz et al. (2000), Freed and Lin (2001) and Masterlark and Wang (2002). With a rigidity μ=30 GPa, the corresponding Maxwell relaxation time is $\mathit{\tau }=\mathit{\eta }/\mathit{\mu }\simeq \mathrm{5}\mathrm{a}$.

Figure 6 Geometry of the Landers (LAN) and Hector Mine (HM), California, faults that generated the 1992 and 1999 earthquakes, respectively. The stars indicate the hypocentres of the seismic events. The labels 1 and 2 identify the asperities on the Landers fault.

We model the 1992 earthquake as a two-mode event 01-10 starting from mode 00. Accordingly, the orbit of the system during mode 00 lies inside the subset H₂ of the sticking region and the state P₁ at the beginning of the earthquake belongs to segment s₂ (Fig. 4). The coordinates of P₁ are

$$\begin{array}{}\text{(59)}& X_{\mathrm{1}}=\mathit{\alpha }{Z}_{\mathrm{1}}+\mathrm{1}-\mathit{\alpha }\mathit{\beta }\mathit{\kappa }U,Y_{\mathrm{1}}=\mathit{\beta }-\mathit{\alpha }{Z}_{\mathrm{1}},Z_{\mathrm{1}},\end{array}$$

with

$$\begin{array}{}\text{(60)}& Z_a\le Z_{\mathrm{1}}\le Z_b,\end{array}$$

where the extreme values Z_a and Z_b correspond to the end points of s₂:

$$\begin{array}{}\text{(61)}& Z_a={\displaystyle \frac{\mathit{\kappa }U(\mathit{\alpha }\mathit{\beta }+\mathrm{1})-\mathrm{1}}{\mathit{\alpha }}},Z_b={\displaystyle \frac{\mathit{\beta }(\mathrm{1}-\mathit{\kappa }U)}{\mathit{\alpha }}}.\end{array}$$

At the end of mode 01, the system is at point P₂ with coordinates

$$\begin{array}{}\text{(62)}& X_{\mathrm{2}}=X_{\mathrm{1}},Y_{\mathrm{2}}=Y_{\mathrm{1}}-\mathit{\beta }\mathit{\kappa }U,Z_{\mathrm{2}}=Z_{\mathrm{1}}-\mathit{\beta }\mathit{\kappa }U,\end{array}$$

where mode 10 starts. As Z₁ varies in the interval given in Eq. (37), an infinite number of points P₂ describe a segment r₂ on the subset Q₁ of the face AECD and parallel to the edge CD. Mode 10 terminates at point P₃ with coordinates

$$\begin{array}{}\text{(63)}& X_{\mathrm{3}}=X_{\mathrm{2}}-\mathit{\kappa }U,Y_{\mathrm{3}}=Y_{\mathrm{2}},Z_{\mathrm{3}}=Z_{\mathrm{2}}+\mathit{\kappa }U.\end{array}$$

Again, as Z₁ varies in the interval given in Eq. (37), there is an infinite number of points P₃ defining another segment q₂ parallel to the edge CD. This segment is situated within the sticking region and crosses the surface Σ for Z₁=Z_c, with $Z_a<Z_c<Z_b$.

Dragoni and Tallarico (2016) studied the 1992 Landers earthquake under the hypothesis of purely elastic coupling between the asperities. Following the authors, we take α=0.1, β=0.5, γ=1.5 and ϵ=0.7, a set of values yielding modelled moment rate and seismic spectrum comparable with the observations. Thus, we have U≃0.546 and κ≃0.52. As for viscoelastic relaxation, it can be characterized by the product VΘ (Amendola and Dragoni, 2013), which can be estimated as

$$\begin{array}{}\text{(64)}& V\mathrm{\Theta }={\displaystyle \frac{\mathit{\kappa }Uv\mathit{\tau }}{u_{\mathrm{1}}}},\end{array}$$

where v=3 cm a⁻¹ is the relative plate velocity at Landers (Wallace, 1990). Accordingly, we have VΘ≃0.007.

Every state P₁ on segment s₂, where the 1992 earthquake began, corresponds to a specific state P₃ on segment q₂, where the 1992 earthquake ended. Exploiting Eq. (39), we can express the coordinates (40) of P₃ as a function of Z₁. Since q₂ crosses the surface Σ, the state P₃ can belong to H₁,H₂ or Σ, in correspondence to $Z_c<Z_{\mathrm{1}}\le Z_b$, $Z_a\le Z_{\mathrm{1}}<Z_c$ and Z₁=Z_c, respectively. In the first case, the next event will start with the failure of asperity 1; in the second case, with the failure of asperity 2; in the third case, with the simultaneous failure of the asperities. With the values of α, β, κ and U listed above, we find $Z_a\simeq -\mathrm{7.02}$, Z_b≃3.58 and Z_c≃0.78. Accordingly, only about one-quarter of segment q₂ lies inside the subset H₁ of the sticking region. Without any further discussion and neglecting the stress perturbation caused by the Hector Mine earthquake, we would infer that future events on the 1992 fault are more likely to start with the failure of asperity 2.

5.1 Stress perturbation by the 1999 Hector Mine earthquake

The 1999 Hector Mine earthquake was generated by right-lateral strike-slip faulting located at 34.59^∘ N, 116.27^∘ W, about 20 km north-east from the Landers fault (Jónsson et al., 2002; Salichon et al., 2004). We characterize the event averaging the data available in the SRCMOD database and assume strike, dip and rake angles of 330, 80 and 180^∘, respectively; a depth of 10 km; and a seismic moment of 6.62×10¹⁹ Nm.

The stress transferred to the asperities at Landers can be evaluated employing the model of Appendix Appendix A, taking

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{1}}=\mathrm{345}{}^{\circ },{\mathit{\varphi }}_{\mathrm{2}}=\mathrm{330}{}^{\circ },{\mathit{\psi }}_{\mathrm{1}}=\mathrm{85}{}^{\circ },{\mathit{\psi }}_{\mathrm{2}}=\mathrm{80}{}^{\circ },\\ \text{(65)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{1}}={\mathit{\lambda }}_{\mathrm{2}}=\mathrm{180}{}^{\circ }.\end{array}$$

As a result, the normal and tangential components of the perturbing stress on asperity 1 are

$$\begin{array}{}\text{(66)}& {\mathit{\sigma }}_{\mathrm{1}n}\simeq \mathrm{0.14}\mathrm{MPa},{\mathit{\sigma }}_{\mathrm{1}t}\simeq \mathrm{0.39}\mathrm{MPa}.\end{array}$$

Accordingly, the static friction on asperity 1 is reduced and right-lateral slip is favoured. As for asperity 2, the components of the perturbing stress are

$$\begin{array}{}\text{(67)}& {\mathit{\sigma }}_{\mathrm{2}n}\simeq \mathrm{0.18}\mathrm{MPa},{\mathit{\sigma }}_{\mathrm{2}t}\simeq -\mathrm{0.17}\mathrm{MPa},\end{array}$$

suggesting that static friction on asperity 2 is reduced and right-lateral slip is inhibited.

We now introduce the effect of the perturbation in the framework of the discrete model. The changes in the tangential forces (2) on the asperities are

$$\begin{array}{}\text{(68)}& \mathrm{\Delta }{F}_{\mathrm{1}}=-{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}t}}{f_{\mathrm{s}\mathrm{1}}}}A,\mathrm{\Delta }{F}_{\mathrm{2}}=-{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{2}t}}{f_{\mathrm{s}\mathrm{1}}}}A.\end{array}$$

The static friction f_s1 on asperity 1 can be evaluated as (Dragoni and Santini, 2012)

$$\begin{array}{}\text{(69)}& f_{\mathrm{s}\mathrm{1}}={\displaystyle \frac{K{u}_{\mathrm{1}}}{\mathit{\kappa }U}},\end{array}$$

where the constant

$$\begin{array}{}\text{(70)}& K={\displaystyle \frac{\mathit{\mu }A}{d}}\end{array}$$

is an expression of the coupling between the asperities and the tectonic plates. With d=80 km (Masterlark and Wang, 2002), it results in $f_{\mathrm{s}\mathrm{1}}/A\simeq \mathrm{7.9}$ MPa. Hence, we have

$$\begin{array}{}\text{(71)}& \mathrm{\Delta }{F}_{\mathrm{1}}\simeq -\mathrm{0.05},\mathrm{\Delta }{F}_{\mathrm{2}}\simeq \mathrm{0.02}.\end{array}$$

From Eqs. (19)–(19), the components of the perturbation vector ΔR are

$$\begin{array}{}\text{(72)}& \mathrm{\Delta }X\simeq \mathrm{0.043},\mathrm{\Delta }Y\simeq -\mathrm{0.016},\mathrm{\Delta }Z\simeq -\mathrm{0.059}.\end{array}$$

As a result, the orientation of ΔR in the state space is characterized by angles $\mathit{\theta }\simeq -\mathrm{0.35}$ rad and $\mathit{\delta }\simeq -\mathrm{0.91}$ rad. The changes in static frictions (23) can be calculated as

$$\begin{array}{}\text{(73)}& \mathrm{\Delta }{\mathit{\beta }}_{\mathrm{1}}=-{\displaystyle \frac{k_{\mathrm{s}}{\mathit{\sigma }}_{\mathrm{1}n}}{f_{\mathrm{s}\mathrm{1}}}}A,\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{2}}=-{\displaystyle \frac{k_{\mathrm{s}}{\mathit{\sigma }}_{\mathrm{2}n}}{f_{\mathrm{s}\mathrm{1}}}}A,\end{array}$$

where k_s is the effective static friction coefficient on asperity 1. Assuming k_s=0.4, we get

$$\begin{array}{}\text{(74)}& \mathrm{\Delta }{\mathit{\beta }}_{\mathrm{1}}\simeq -\mathrm{0.0073},\mathrm{\Delta }{\mathit{\beta }}_{\mathrm{2}}\simeq -\mathrm{0.0092}.\end{array}$$

Finally, from Eqs. (29) and (29), the changes in Coulomb forces on the asperities are

$$\begin{array}{}\text{(75)}& \mathrm{\Delta }{F}_{\mathrm{1}}^{\mathrm{C}}\simeq \mathrm{0.057},\mathrm{\Delta }{F}_{\mathrm{2}}^{\mathrm{C}}\simeq -\mathrm{0.012}.\end{array}$$

At the time of the Hector Mine earthquake, the Landers fault was at a state P₄ resulting from the post-seismic evolution of any of the possible states P₃∈q₂ where the 1992 event ended. The coordinates of P₄ can be calculated from the solution to the equations of mode 00 given by Dragoni and Lorenzano (2015) and taking into account that the time interval $\stackrel{\mathrm{̃}}{t}$ elapsed between the 1992 Landers and 1999 Hector Mine earthquakes amounts to about 7.3 years:

$$\begin{array}{}\text{(76)}& X_{\mathrm{4}}=X_{\mathrm{3}}+V\mathrm{\Theta }\stackrel{\mathrm{̃}}{T},Y_{\mathrm{4}}=Y_{\mathrm{3}}+V\mathrm{\Theta }\stackrel{\mathrm{̃}}{T},Z_{\mathrm{4}}=Z_{\mathrm{3}}{e}^{-\stackrel{\mathrm{̃}}{T}},\end{array}$$

where

$$\begin{array}{}\text{(77)}& \stackrel{\mathrm{̃}}{T}={\displaystyle \frac{\stackrel{\mathrm{̃}}{t}}{\mathit{\tau }}}\approx \mathrm{1.5}.\end{array}$$

Making use of Eqs. (39) and (40), we can express the coordinates of P₄ as a function of $Z_{\mathrm{1}}\in [Z_a,Z_b]$. Accordingly, there is an infinite number of points P₄ defining a vector t₂ inside the sticking region. At $T=\stackrel{\mathrm{̃}}{T}$, the perturbation vector ΔR moves every state P₄ to a new state $P_{\mathrm{4}}^{\prime }$ with coordinates

$$\begin{array}{}\text{(78)}& X_{\mathrm{4}}^{\prime }=X_{\mathrm{4}}+\mathrm{\Delta }X,Y_{\mathrm{4}}^{\prime }=Y_{\mathrm{4}}+\mathrm{\Delta }Y,Z_{\mathrm{4}}^{\prime }=Z_{\mathrm{4}}+\mathrm{\Delta }Z,\end{array}$$

which can be expressed as a function of $Z_{\mathrm{1}}\in [Z_a,Z_b]$. As a result, a new vector ${\mathit{t}}_{\mathrm{2}}^{\prime }$ identifies the state of the Landers fault after the Hector Mine earthquake.

In order to characterize the effect of the perturbation, let us consider the difference ΔF^C defined in Eq. (29): from Eq. (51), we get $\mathrm{\Delta }{F}^{\mathrm{C}}\simeq -\mathrm{0.069}$. Since ΔF^C<0, we conclude that the stress perturbation is such that states P₄∈H₁ are moved to ${\mathbf{H}}_{\mathrm{1}}^{\prime }$, the state P₄∈Σ enters ${\mathbf{H}}_{\mathrm{1}}^{\prime }$, states P₄∈H₂ are shifted towards the Σ surface and some of them enter ${\mathbf{H}}_{\mathrm{1}}^{\prime }$. Specifically, we find that $P_{\mathrm{4}}^{\prime }$ belongs to ${\mathbf{H}}_{\mathrm{1}}^{\prime }$, ${\mathbf{H}}_{\mathrm{2}}^{\prime }$ and Σ^′ in correspondence to $Z_c^{\prime }<Z_{\mathrm{1}}\le Z_b$, $Z_a\le Z_{\mathrm{1}}<Z_c^{\prime }$ and $Z_{\mathrm{1}}=Z_c^{\prime }$, with $Z_c^{\prime }\simeq \mathrm{0.50}$. On the whole, we can draw the preliminary conclusion that the stress perturbation is such that future events on the Landers fault starting with the slip of asperity 1 are favoured over events starting with the slip of asperity 2. A deeper discussion is provided in the following.

5.2 Constraints due to the seismic history to date

In order to improve our knowledge on the state that gave rise to the 1992 Landers earthquake and on the possible future events generated by that fault, we exploit the seismic history between 1999 and the present date. After the perturbation caused by the Hector Mine earthquake, the inter-seismic time $T_{\mathrm{is}}^{\prime }$ of the Landers fault can be calculated from Eq. (30) and Eq. (33) for states $P_{\mathrm{4}}^{\prime }$ belonging to ${\mathbf{H}}_{\mathrm{1}}^{\prime }$ and ${\mathbf{H}}_{\mathrm{2}}^{\prime }$, respectively:

$$\begin{array}{}\text{(79)}& T_{\mathrm{is}}^{\prime }=\left\{\begin{array}{l}\mathrm{\Theta }W\left({\mathit{\gamma }}_{\mathrm{1}}^{\prime }\right)+{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{1}}-X_{\mathrm{4}}^{\prime }}{V}},Z_c^{\prime }<Z_{\mathrm{1}}\le Z_b\\ \mathrm{\Theta }W\left({\mathit{\gamma }}_{\mathrm{2}}^{\prime }\right)+{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{2}}-Y_{\mathrm{4}}^{\prime }}{V}},Z_a\le Z_{\mathrm{1}}<Z_c^{\prime }\end{array}\right.,\end{array}$$

where

$$\begin{array}{}\text{(80)}& {\mathit{\gamma }}_{\mathrm{1}}^{\prime }={\displaystyle \frac{\mathit{\alpha }{Z}_{\mathrm{4}}^{\prime }}{V\mathrm{\Theta }}}{e}^{-\frac{{\mathit{\beta }}_{\mathrm{1}}-X_{\mathrm{4}}^{\prime }}{V\mathrm{\Theta }}},{\mathit{\gamma }}_{\mathrm{2}}^{\prime }=-{\displaystyle \frac{\mathit{\alpha }{Z}_{\mathrm{4}}^{\prime }}{V\mathrm{\Theta }}}{e}^{-\frac{{\mathit{\beta }}_{\mathrm{2}}-Y_{\mathrm{4}}^{\prime }}{V\mathrm{\Theta }}}.\end{array}$$

Since no earthquakes have been produced by the Landers fault after the occurrence of the Hector Mine event, up to year 2016, we can exclude the states on the segment s₂ yielding an expected inter-seismic time (Eq. 55) shorter than or equal to $t_{\mathrm{is}}^{\prime }=$ 17 years. The requirement

$$\begin{array}{}\text{(81)}& T_{\mathrm{is}}^{\prime }>{\displaystyle \frac{t_{\mathrm{is}}^{\prime }}{\mathit{\tau }}}\mathrm{\Theta }\approx \mathrm{3.5}\mathrm{\Theta }\end{array}$$

is satisfied by states on segment s₂ in the subset ${\stackrel{\mathrm{̃}}{Z}}_a\le Z_{\mathrm{1}}\le {\stackrel{\mathrm{̃}}{Z}}_b$, with ${\stackrel{\mathrm{̃}}{Z}}_a\simeq -\mathrm{1.17}$ and ${\stackrel{\mathrm{̃}}{Z}}_b\simeq \mathrm{2.19}$.

As a consequence, we can constrain the admissible states on the segment t₂. A comparison between the intervals $[{\stackrel{\mathrm{̃}}{Z}}_a,Z_c]$ and $[Z_c,{\stackrel{\mathrm{̃}}{Z}}_b]$ points out that more than one-half of the acceptable subset of t₂ belongs to H₂. Hence, before the stress perturbation caused by the Hector Mine earthquake, future events on the 1992 Landers fault were more likely to start with the failure of asperity 2.
In turn, the refinement of t₂ limits the acceptable states on the segment ${\mathit{t}}_{\mathrm{2}}^{\prime }$. From the amplitude of the intervals $[{\stackrel{\mathrm{̃}}{Z}}_a,Z_c^{\prime }]$ and $[Z_c^{\prime },{\stackrel{\mathrm{̃}}{Z}}_b]$, we deduce that the acceptable subset of ${\mathit{t}}_{\mathrm{2}}^{\prime }$ is almost equally divided between ${\mathbf{H}}_{\mathrm{1}}^{\prime }$ and ${\mathbf{H}}_{\mathrm{2}}^{\prime }$. Therefore, if we consider the influence of the Hector Mine earthquake on future events generated by the 1992 Landers fault, we conclude that the stress perturbation yielded homogenization in the probability of events starting with the failure of asperity 1 or asperity 2. This result is in agreement with the observation that the perturbation vector ΔR shifted the whole segment t₂ towards the subset ${\mathbf{H}}_{\mathrm{1}}^{\prime }$ of the sticking region.

These conclusions would have to be reconsidered if new stress perturbations from neighbouring faults were to affect the post-seismic evolution of the Landers fault in the future. In addition, if no earthquakes were to be observed for some time on the Landers fault, the refining procedure discussed above could be repeated and the admissible subsets of segments s₂, t₂ and ${\mathit{t}}_{\mathrm{2}}^{\prime }$ could be constrained with further precision.

5.3 Effects of the stress perturbation on future earthquakes

Finally, we discuss the features of the next seismic event generated by the 1992 Landers fault, highlighting the changes due to the Hector Mine earthquake.

Every state P₁∈s₂ where the 1992 earthquake began corresponds to a particular state P₄∈t₂ and $P_{\mathrm{4}}^{\prime }\in {\mathit{t}}_{\mathrm{2}}^{\prime }$ before and after the stress perturbation associated with the Hector Mine earthquake, respectively. Since the segment t₂ intersects the surface Σ, the state P₄ can belong to H₁, H₂ or Σ (Fig. 3), thus affecting the asperity that will fail the first at the beginning of the next earthquake on the fault. In the first case, the next event will start with the failure of asperity 1, in the second case with the failure of asperity 2, in the third case with the simultaneous failure of the asperities. Analogous considerations hold for states $P_{\mathrm{4}}^{\prime }$ in ${\mathbf{H}}_{\mathrm{1}}^{\prime },{\mathbf{H}}_{\mathrm{2}}^{\prime }$ and Σ^′, respectively.

The number and the sequence of dynamic modes in the earthquake depend on the sub-interval of Z₁ considered. The details are summarized in Table 3 for both the unperturbed and perturbed cases. Taking these specifics into account and referring to Tables 1 and 2, we evaluate the seismic moments M₀ and $M_{\mathrm{0}}^{\prime }$ associated with the expected future earthquake on the 1992 fault before and after the Hector Mine earthquake, respectively. In Fig. 7, we show the difference

$$\begin{array}{}\text{(82)}& \mathrm{\Delta }{M}_{\mathrm{0}}=M_{\mathrm{0}}^{\prime }-M_{\mathrm{0}}\end{array}$$

as a function of $Z_{\mathrm{1}}\in [{\stackrel{\mathrm{̃}}{Z}}_a,{\stackrel{\mathrm{̃}}{Z}}_b]$. Owing to the translation imposed to the segment t₂ by the perturbation vector ΔR, the sign of ΔM₀ changes across the different sub-intervals of Z₁. The energy released by the earthquake is increased for $Z_{\mathrm{1}}\in [\mathrm{0.43},\mathrm{0.71}]$, while it is reduced elsewhere.

Another significant result of the stress perturbation concerns the variation in the inter-seismic time before the next seismic event. As in Sect. 5.2, we consider the post-seismic evolution from 1999 onwards and set the origin of times at the occurrence of the Hector Mine earthquake. The expected inter-seismic time T_is prior to the stress perturbation can be calculated from Eqs. (29) and (32) for states P₄ belonging to H₁ and H₂, respectively:

$$\begin{array}{}\text{(83)}& T_{\mathrm{is}}=\left\{\begin{array}{l}\mathrm{\Theta }W\left({\mathit{\gamma }}_{\mathrm{1}}\right)+{\displaystyle \frac{\mathrm{1}-X_{\mathrm{4}}}{V}},Z_c<Z_{\mathrm{1}}\le {\stackrel{\mathrm{̃}}{Z}}_b\\ \mathrm{\Theta }W\left({\mathit{\gamma }}_{\mathrm{2}}\right)+{\displaystyle \frac{\mathit{\beta }-Y_{\mathrm{4}}}{V}},{\stackrel{\mathrm{̃}}{Z}}_a\le Z_{\mathrm{1}}<Z_c\end{array}\right.,\end{array}$$

where

$$\begin{array}{}\text{(84)}& {\mathit{\gamma }}_{\mathrm{1}}={\displaystyle \frac{\mathit{\alpha }{Z}_{\mathrm{4}}}{V\mathrm{\Theta }}}{e}^{-\frac{\mathrm{1}-X_{\mathrm{4}}}{V\mathrm{\Theta }}},{\mathit{\gamma }}_{\mathrm{2}}=-{\displaystyle \frac{\mathit{\alpha }{Z}_{\mathrm{4}}}{V\mathrm{\Theta }}}{e}^{-\frac{\mathit{\beta }-Y_{\mathrm{4}}}{V\mathrm{\Theta }}}.\end{array}$$

The inter-seismic time $T_{\mathrm{is}}^{\prime }$ after the stress perturbation has been given in Eq. (55). The difference

$$\begin{array}{}\text{(85)}& \mathrm{\Delta }T=T_{\mathrm{is}}^{\prime }-T_{\mathrm{is}}\end{array}$$

is shown in Fig. 8 as a function of $Z_{\mathrm{1}}\in [{\stackrel{\mathrm{̃}}{Z}}_a,{\stackrel{\mathrm{̃}}{Z}}_b]$. For states P₄∈H₁ corresponding to $P_{\mathrm{4}}^{\prime }\in {\mathbf{H}}_{\mathrm{1}}^{\prime }$ and states P₄∈H₂ corresponding to $P_{\mathrm{4}}^{\prime }\in {\mathbf{H}}_{\mathrm{2}}^{\prime }$, this difference coincides with Eqs. (31) and (34), respectively.

Table 3 Future earthquakes generated by the 1992 Landers, California, fault, as functions of the variable Z₁ describing the initial state of the 1992 event, with $Z_{\mathrm{1}}\in [{\stackrel{\mathrm{̃}}{Z}}_a,{\stackrel{\mathrm{̃}}{Z}}_b]=[-\mathrm{1.16},\mathrm{2.19}]$. The results predicted by the model before and after the stress perturbation associated with the 1999 Hector Mine, California, earthquake are shown. The values $Z_{\mathrm{1}}=Z_c=\mathrm{0.78}$ and $Z_{\mathrm{1}}=Z_c^{\prime }=\mathrm{0.50}$ correspond to the largest possible earthquakes before and after the stress perturbation, respectively.

Figure 7 Change in the seismic moment released during the next event on the 1992 Landers, California, fault, as a result of the stress perturbation due to the 1999 Hector Mine, California, earthquake. On the horizontal axis, the variable Z₁ describing the initial state of the 1992 event. The kinds of seismic event predicted by the model, corresponding to different intervals of Z₁, are listed in Table 3. The values Z₁=Z_c and $Z_{\mathrm{1}}=Z_c^{\prime }$ correspond to the largest possible earthquakes before and after the stress perturbation, respectively.

Some peculiar features stand out. First, we notice that, for all states P₄∈H₂ corresponding to $P_{\mathrm{4}}^{\prime }\in {\mathbf{H}}_{\mathrm{2}}^{\prime }$, that is, for $Z_{\mathrm{1}}\in [{\stackrel{\mathrm{̃}}{Z}}_a,Z_c^{\prime }]$, the inter-seismic time is increased by the stress perturbation, in agreement with the inhibiting effect on asperity 2 suggested by Eq. (51). However, Eq. (51) suggests that the failure of asperity 1 is promoted, but this is not verified by all states $P_{\mathrm{4}}^{\prime }\in {\mathbf{H}}_{\mathrm{1}}^{\prime }$, that is, for $Z_{\mathrm{1}}\in [Z_c^{\prime },{\stackrel{\mathrm{̃}}{Z}}_b]$. In fact, the inter-seismic time is reduced only for $Z_{\mathrm{1}}\in (\mathrm{0.53},{\stackrel{\mathrm{̃}}{Z}}_b]$, while it is increased for $Z_{\mathrm{1}}\in [Z_c^{\prime },\mathrm{0.53})$. In the particular case Z₁=0.53, there is no change in the inter-seismic time. This is a remarkable result, showing that the presence of viscoelastic relaxation at the time of the stress perturbation entails the unpredictability of the consequent influence in terms of anticipation or delay of future earthquakes, on the basis of the sole knowledge of the change in Coulomb stress.

Figure 8 Change in the inter-seismic time before the next event on the 1992 Landers, California, fault, as a result of the stress perturbation due to the 1999 Hector Mine, California, earthquake. On the horizontal axis, the variable Z₁ describing the initial state of the 1992 event. The kinds of seismic event predicted by the model, corresponding to different intervals of Z₁, are listed in Table 3. The values Z₁=Z_c and $Z_{\mathrm{1}}=Z_c^{\prime }$ correspond to the largest possible earthquakes before and after the stress perturbation, respectively.

At the occurrence of the next earthquake produced by the Landers fault, the number and sequence of dynamic modes involved and the energy released will reveal more about the state of the system, thus allowing a further refinement of the specific conditions that gave rise to the 1992 event.

6 Conclusions

We considered a plane fault embedded in a shear zone, subject to a uniform strain rate owing to tectonic loading. The fault is characterized by the presence of two asperities with equal areas and different frictional resistance. The coseismic static stress field due to earthquakes produced by the fault is relaxed by viscoelastic deformation in the asthenosphere.

The fault was treated as a discrete dynamical system with three degrees of freedom: the slip deficits of the asperities and the variation of their difference due to viscoelastic deformation. The dynamics of the system was described in terms of one sticking mode and three slipping modes. In the sticking mode, the orbit of the system lies in a convex hexahedron in the space of the state variables, while the number and the sequence of slipping modes during a seismic event are determined by the particular state of the system at the beginning of the inter-seismic interval preceding the event. The amount of slip of the asperities and the energy released during an earthquake generated by the fault can be predicted accordingly.

The effect of stress transfer due to earthquakes on neighbouring faults was studied in terms of a perturbation vector yielding changes to the state of the system, its sticking region and the energy released during a subsequent seismic event. The specific effect on the evolution of the fault is related with the orientation of this vector in the state space.

We investigated the interplay between the ongoing viscoelastic relaxation on the fault and a stress perturbation imposed during an inter-seismic interval. Following a stress perturbation due to earthquakes on neighbouring faults, an increase in the Coulomb stress associated with a given asperity directly yields the anticipation of the slip of that asperity, and vice versa, if a purely elastic rheology is assumed for the receiving fault (Dragoni and Piombo, 2015). According to the present model, this property no longer holds if the change in Coulomb stress occurs while viscoelastic relaxation is taking place on the receiving fault. In fact, even if the change in the inter-seismic intervals of the asperities can still be evaluated from a theoretical point of view, the specific effect of the stress perturbation could be univocally inferred only if the particular states of the fault at the time of the stress perturbation and right after it were known. The information on the change in Coulomb stress on the fault do not suffice any more.

We applied the model to the stress perturbation imposed by the 1999 Hector Mine, California, earthquake to the fault that caused the 1992 Landers, California, earthquake, which was due to the failure of two asperities and was followed by significant viscoelastic relaxation. We modelled the 1992 Landers earthquake as a two-mode event associated with the separate slip of the asperities and showed how the event is compatible with a number of possible initial states of the fault, which can be screened on the basis of the seismic history to date. The details of the stress transfer associated with the 1999 Hector Mine earthquake were calculated using the relative positions and faulting styles of the two faults as a starting point. We discussed the effect of the stress perturbation, pointing out the complexity of its influence on the possible future events generated by the 1992 Landers fault in terms of the associated energy release, the sequence of dynamic modes involved and the duration of the inter-seismic interval. Specifically, we showed that the consequences of the 1999 Hector Mine earthquake on the post-seismic evolution of the 1992 Landers fault depend on the specific state of the Landers fault at the time of the 1999 earthquake and immediately after it, even if the variations in the Coulomb stress on the asperities at Landers are known. On the whole, the application allowed the exemplification of the critical unpredictability of the effect of a stress perturbation occurring while viscoelastic relaxation is taking place.

Another source of complication may be represented by the interaction between viscoelastic relaxation and stable creep on the fault. This problem is beyond the scope of the present work, but it may be object of future research by combining elements of the present model with the one of Dragoni and Lorenzano (2017).

Appendix A: Estimate of the stress perturbation

We consider two plane faults, namely fault 1 and fault 2, embedded in an infinite, homogeneous and isotropic Poisson medium of rigidity μ (Fig. A1). Following the slip of fault 1 (perturbing fault), stress is transferred to fault 2 (receiving fault). We calculate the normal traction σ_n and the tangential traction in the direction of slip σ_t transferred to the receiving fault, estimated as the average value at its centre.

We define a coordinate system $(x,y,z)$ with axes corresponding with the directions of dip, strike and normal on fault 1, respectively. Fault 1 lies on the plane z=0 and its centre is in the origin of the coordinate system. Accordingly, the unit vector perpendicular to fault 1 is $n_{\mathrm{1}i}=(\mathrm{0},\mathrm{0},\mathrm{1})$. We call ϕ₁, ψ₁ and λ₁ the strike, dip and rake angles of fault 1, respectively. The slip direction of fault 1 is then given by

$$\begin{array}{}\text{(A1)}& m_{\mathrm{1}i}=\left(-\sin {\mathit{\lambda }}_{\mathrm{1}},\cos {\mathit{\lambda }}_{\mathrm{1}},\mathrm{0}\right).\end{array}$$

Fault 2 is characterized by strike and dip angles ϕ₂ and ψ₂, respectively. Accordingly, the unit vector perpendicular to fault 2 is given by

$$\begin{array}{}\text{(A2)}& n_{\mathrm{2}i}=\left(\sin \mathrm{\Delta }\mathit{\psi }\cos \mathrm{\Delta }\mathit{\varphi },-\sin \mathrm{\Delta }\mathit{\psi }\sin \mathrm{\Delta }\mathit{\varphi },\cos \mathrm{\Delta }\mathit{\psi }\right),\end{array}$$

where

$$\begin{array}{}\text{(A3)}& \mathrm{\Delta }\mathit{\varphi }={\mathit{\varphi }}_{\mathrm{2}}-{\mathit{\varphi }}_{\mathrm{1}},\mathrm{\Delta }\mathit{\psi }={\mathit{\psi }}_{\mathrm{2}}-{\mathit{\psi }}_{\mathrm{1}}.\end{array}$$

Let λ₂ be the preferred rake angle on fault 2, correlated with the orientation of tectonic loading: λ₂=0^∘ for left-lateral strike slip, λ₂=180^∘ for right-lateral strike slip, ${\mathit{\lambda }}_{\mathrm{2}}=-\mathrm{90}$^∘ for normal dip slip and λ₂=90^∘ for reverse dip slip. The corresponding slip direction is

$$\begin{array}{}\text{(A4)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{2}x}=\cos {\mathit{\lambda }}_{\mathrm{2}}\sin \mathrm{\Delta }\mathit{\varphi }-\sin {\mathit{\lambda }}_{\mathrm{2}}\cos \mathrm{\Delta }\mathit{\psi }\cos \mathrm{\Delta }\mathit{\varphi },\text{(A5)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{2}y}=\cos {\mathit{\lambda }}_{\mathrm{2}}\cos \mathrm{\Delta }\mathit{\varphi }+\sin {\mathit{\lambda }}_{\mathrm{2}}\cos \mathrm{\Delta }\mathit{\psi }\sin \mathrm{\Delta }\mathit{\varphi },\text{(A6)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{2}z}=\sin {\mathit{\lambda }}_{\mathrm{2}}\sin \mathrm{\Delta }\mathit{\psi }.\end{array}$$

We name (E_i,N_i) and D_i the UTM (Universal Transverse Mercator) coordinates and depths of the centres of the faults, respectively. In our reference system, the coordinates of the centre of fault 2 are identified by the following three steps:

1. placing the origin at the centre of fault 1:

   $$\begin{array}{}\text{(A7)}& x^{\prime }=E_{\mathrm{2}}-E_{\mathrm{1}},y^{\prime }=N_{\mathrm{2}}-N_{\mathrm{1}},z^{\prime }=D_{\mathrm{2}}-D_{\mathrm{1}}\end{array}$$

2. clockwise rotation about the z axis by the angle ϕ₁:

   $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{x}^{\prime \prime }=x^{\prime }\cos {\mathit{\varphi }}_{\mathrm{1}}-y^{\prime }\sin {\mathit{\varphi }}_{\mathrm{1}}{y}^{\prime \prime }=x^{\prime }\sin {\mathit{\varphi }}_{\mathrm{1}}+y^{\prime }\cos {\mathit{\varphi }}_{\mathrm{1}},\\ \text{(A8)}& {\displaystyle }& {\displaystyle }{z}^{\prime \prime }=z^{\prime }\end{array}$$

3. counterclockwise rotation about the y axis by the angle ψ₁:

   $$\begin{array}{ll}{\displaystyle }& {\displaystyle }x=x^{\prime \prime }\cos {\mathit{\psi }}_{\mathrm{1}}-z^{\prime \prime }\sin {\mathit{\psi }}_{\mathrm{1}},y=y^{\prime \prime },\\ \text{(A9)}& {\displaystyle }& {\displaystyle }z=x^{\prime \prime }\sin {\mathit{\psi }}_{\mathrm{1}}+z^{\prime \prime }\cos {\mathit{\psi }}_{\mathrm{1}}.\end{array}$$

Figure A1 Geometry of the model employed to study the stress transfer between neighbouring faults. Fault 1 is the perturbing fault, while fault 2 is the receiving fault. The coordinates $(E,N,D)$ are the UTM coordinates and depth of the centres of the faults, respectively, whereas the axes $(x,y,z)$ correspond with the directions of dip, strike and normal on fault 1, respectively. The angles ϕ and ψ are the strike and dip angles of the faults, respectively.

The perturbing fault is treated as a point-like dislocation source (a double couple of forces) located at the origin. This is a good approximation for non-overlapping regions (Dragoni and Lorenzano, 2016). Let m₀ be the scalar seismic moment of the dislocation. The ith component of the static displacement field generated by the slip of fault 1 is

$$\begin{array}{}\text{(A10)}& u_i=-M_{jk}{G}_{ij,k},\end{array}$$

where M_ij is the moment tensor associated with the dislocation source

$$\begin{array}{}\text{(A11)}& M_{ij}=m_{\mathrm{0}}(m_{\mathrm{1}i}{n}_{\mathrm{1}j}+m_{\mathrm{1}j}{n}_{\mathrm{1}i})\end{array}$$

and G_ij is the Somigliana tensor

$$\begin{array}{}\text{(A12)}& G_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{8}\mathit{\pi }\mathit{\mu }}}\left({\displaystyle \frac{\mathrm{2}}{r}}{\mathit{\delta }}_{ij}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}r{,}_{ij}\right),\end{array}$$

with

$$\begin{array}{}\text{(A13)}& r=\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}+z^{\mathrm{2}}}.\end{array}$$

The components of the stress field are given by

$$\begin{array}{}\text{(A14)}& {\mathit{\sigma }}_{ij}=\mathit{\mu }(e_{kk}{\mathit{\delta }}_{ij}+\mathrm{2}{e}_{ij}),\end{array}$$

where e_ij is the strain field associated with the displacement field (A5). Finally, the normal traction σ_n and the tangential traction in the direction of slip σ_t on fault 2 are

$$\begin{array}{}\text{(A15)}& {\mathit{\sigma }}_{\mathrm{n}}={\mathit{\sigma }}_{ij}{n}_{\mathrm{2}i}{n}_{\mathrm{2}j},{\mathit{\sigma }}_{\mathrm{t}}={\mathit{\sigma }}_{ij}{m}_{\mathrm{2}i}{n}_{\mathrm{2}j}.\end{array}$$

The signs of σ_n and σ_t define the effect of the stress transfer on fault 2. If σ_n>0, the amount of compressional stress on the receiving fault is reduced, and vice versa. If σ_t>0, the slip of the receiving fault is promoted, and vice versa.