Stress states and moment rates of a two-asperity fault in the presence of viscoelastic relaxation

Abstract

A fault containing two asperities with different strengths is considered. The fault is embedded in a shear zone subject to a constant strain rate by the motions of adjacent tectonic plates. The fault is modelled as a discrete dynamical system where the average values of stress, friction and slip on each asperity are considered. The state of the fault is described by three variables: the slip deficits of the asperities and the viscoelastic deformation. The system has four dynamic modes, for which analytical solutions are calculated. The relationship between the state of the fault before a seismic event and the sequence of slipping modes in the event is enlightened. Since the moment rate depends on the number and sequence of slipping modes, the knowledge of the source function of an earthquake constrains the orbit of the system in the phase space. If the source functions of a larger number of consecutive earthquakes were known, the orbit could be constrained more and more and its evolution could be predicted with a smaller uncertainty. The model is applied to the 1964 Alaska earthquake, which was the effect of the failure of two asperities and for which a remarkable post-seismic relaxation has been observed in the subsequent decades. The evolution of the system after the 1964 event depends on the state from which the event was originated, that is constrained by the observed moment rate. The possible durations of the interseismic interval and the possible moment rates of the next earthquake are calculated as functions of the initial state.

Introduction

Many aspects of fault dynamics can be reproduced by asperity models (Lay et al., 1982; Scholz, 1990), assuming that one or more regions of the fault have a much higher friction than the adjacent regions. Several large and medium-size earthquakes that occurred in the last decades were the result of the failure of two distinct asperities, such as the 1964 Alaska earthquake (Christensen and Beck, 1994), the 1995 Kobe earthquake (Kikuchi and Kanamori, 1996), the 2004 Parkfield earthquake (Johanson et al., 2006) and the 2010 Maule, Chile, earthquake (Delouis et al., 2010).

In the framework of an asperity model, the evolution of asperities in terms of stress accumulation, seismic slip and mutual stress transfer plays a key role. Therefore, the dynamical behaviour of a fault can be fruitfully investigated by means of discrete models describing the state of asperities (e.g. Ruff, 1992; Rice, 1993; Turcotte, 1997). An advantage associated with a finite number of degrees of freedom is that we can predict the evolution of the system at long term by calculating its orbit in the phase space.

A discrete fault model with two asperities was originally proposed by Nussbaum and Ruina (1987) and further investigated by Huang and Turcotte (1990, 1992), McCloskey and Bean (1992) and others. Dragoni and Santini (2012, 2014) solved analytically the equations of motion in the case of a two-asperity fault with different strengths in an elastic medium.

In the long-term evolution of a fault, the rheological properties of Earth's lithosphere play an important role. Lithospheric rocks are not perfectly elastic but have a certain degree of anelasticity (Carter, 1976; Kirby and Kronenberg, 1987; Ranalli, 1995; Nishimura and Thatcher, 2003; Bürgmann and Dresen, 2008). As a consequence, the static stress fields produced by fault dislocations undergo a certain amount of relaxation during the interseismic intervals, which alters the stress distribution on faults and modifies the occurrence times of seismic events (Kusznir, 1991; Kenner and Segall, 2000; Smith and Sandwell, 2006; Piombo et al., 2007; Ding and Lin, 2014).

A preliminary study of the effects of viscoelastic relaxation on a fault containing two asperities was made by Amendola and Dragoni (2013), in the case of asperities with the same frictional strength. It was shown that the stresses on the asperities increase non-linearly during the interseismic intervals, although the tectonic loading takes place at constant rate. As a consequence, earthquakes are anticipated or delayed with respect to the case of an elastic medium. In addition, the stress rate is different for the two asperities, so that the stress distribution changes during loading and the asperity subject to the greater stress at a given instant in time is not necessarily the first one to fail in the next earthquake.

The present paper generalizes Amendola and Dragoni (2013) in that it considers two asperities with different strengths and a larger set of possible states for the fault in the interseismic intervals. We investigate which subsets of states drive the system to the failure of one asperity or both. Whether the failure starts at one asperity or the other has consequences on the position of the earthquake focus as well as on its source function and seismic moment.

The model is applied to the 1964 Alaska earthquake, for which a sufficiently long time interval has elapsed to allow for observation of a remarkable post-seismic relaxation (Zweck et al., 2002). The moment rate of the earthquake was modelled by Dragoni and Santini (2012), showing that it can be approximately represented as a two-mode event with the consecutive failure of the two asperities. We study the subsequent evolution of the system in the presence of viscoelastic relaxation and calculate the duration of the interseismic interval and the possible source functions of the next earthquake.

The model

We consider a plane fault with two asperities of equal areas, that we name asperity 1 and 2, respectively (Fig. 1). Following Amendola and Dragoni (2013), all quantities are expressed in nondimensional form. We assume that the fault is embedded in a homogeneous and isotropic shear zone, subject to a uniform strain rate by the motion of two tectonic plates at relative velocity $V$. The rheological properties of the lithosphere are taken into account by assuming that coseismic stresses are relaxed with a characteristic Maxwell time $\mathrm{\Theta }$.

The fault model. The state of the fault is described by the slip deficits $X\left(T\right)$ and $Y\left(T\right)$ of the asperities and by the viscoelastic deformation $Z\left(T\right)$.

[Figure omitted. See PDF]

We do not determine stress, friction and slip at every point of the fault but, instead, the average values of these quantities on each asperity. We define the slip deficit of an asperity at a certain instant $T$ in time as the slip that the asperity should undergo in order to recover the relative plate displacement occurred up to time $T$.

The state of the fault is described by three variables: $X$, $Y$ and $Z$, where $X$ and $Y$ are the slip deficits of asperities 1 and 2, respectively, while $Z$ is viscoelastic deformation. Accordingly, the asperities are subject to tangential forces: $$F_{\mathrm{1}}=-X+\mathit{\alpha }Z,F_{\mathrm{2}}=-Y-\mathit{\alpha }Z,$$ where $\mathit{\alpha }$ is a coupling constant and the terms $\pm \mathit{\alpha }Z$ are the contribution of stress transfer between the asperities in the presence of viscoelastic deformation. The couple ($F_{\mathrm{1}}$, $F_{\mathrm{2}}$) yields the stress state of the fault.

The forces $F_{\mathrm{1}}$ and $F_{\mathrm{2}}$ are defined as the forces that act on the asperities in the slip direction; therefore, they are valid for any source mechanism (strike-slip, dip-slip or other). In Eq. (1), the terms $-X$ and $-Y$ represent the tectonic loading of asperities 1 and 2, respectively, and have the same sign for both asperities. In the expression for $F_{\mathrm{1}}$, the term $\mathit{\alpha }Z$ is the force applied to asperity 1 by the past motions of asperities, in the presence of viscoelastic relaxation. Analogously, in the expression for $F_{\mathrm{2}}$, the term $-\mathit{\alpha }Z$ is the force applied to asperity 2 by the past motions of asperities.

Fault slip is governed by friction, which is best described by the rate-and-state friction laws (Ruina, 1983; Dieterich, 1994). According to the premise, we use a simpler law assuming that the asperities are characterized by constant static frictions and consider the average values of dynamic frictions during fault slip. We assume that the static friction of asperity 2 is a fraction $\mathit{\beta }$ of that of asperity 1 and that dynamic frictions are a fraction $\mathit{ϵ}$ of static frictions.

If we call $f_{\text{s1}}$ and $f_{\text{d1}}$ the static and the dynamic frictions of asperity 1, respectively, and $f_{\text{s2}}$ and $f_{\text{d2}}$ the static and the dynamic frictions of asperity 2, we define $$\mathit{ϵ}=\frac{f_{\text{d1}}}{f_{\text{s1}}}=\frac{f_{\text{d2}}}{f_{\text{s2}}}$$ and $$\mathit{\beta }=\frac{f_{\text{s2}}}{f_{\text{s1}}}=\frac{f_{\text{d2}}}{f_{\text{d1}}}.$$ Hence the system is described by the five parameters $\mathit{\alpha }$, $\mathit{\beta }$, $\mathit{ϵ}$, $\mathrm{\Theta }$ and $V$, with $\mathit{\alpha }$ $>$ 0, 0 $<$ $\mathit{\beta }$ $<$ 1, 0 $<$ $\mathit{ϵ}$ $<$ 1, $\mathrm{\Theta }$ $>$ 0, and $V$ $>$ 0. From these parameters we can define a slip $$U=\mathrm{2}\frac{\mathrm{1}-\mathit{ϵ}}{\mathrm{1}+\mathit{\alpha }}$$ and the frequencies $$\mathit{\omega }=\sqrt{\mathrm{1}+\mathit{\alpha }},\mathrm{\Omega }=\sqrt{\mathrm{1}+\mathrm{2}\mathit{\alpha }}$$ that will appear in the solutions. The system is subject to the additional constraint $$X\ge \mathrm{0},Y\ge \mathrm{0}$$ that excludes overshooting during fault slip. Forces are expressed in terms of the static friction of asperity 1, so that the conditions for the failure of asperities 1 and 2 are, respectively, $$F_{\mathrm{1}}=-\mathrm{1},F_{\mathrm{2}}=-\mathit{\beta }$$ or, from Eq. (), $$\begin{array}{}& & X-\mathit{\alpha }Z-\mathrm{1}=\mathrm{0},& & Y+\mathit{\alpha }Z-\mathit{\beta }=0.\end{array}$$ These are the equations of two planes in the space $XYZ$, that we call ${\mathrm{\Pi }}_{\mathrm{1}}$ and ${\mathrm{\Pi }}_{\mathrm{2}}$, respectively.

The dynamics of the system has four different modes: a sticking mode, corresponding to stationary asperities (mode 00), and three slipping modes, corresponding to the failure of asperity 1 (mode 10), the failure of asperity 2 (mode 01), and the simultaneous failure of both asperities (mode 11). Each mode is described by a different system of differential equations.

In mode 00, the velocities $\dot{X}$, $\dot{Y}$ and $\dot{Z}$ are negligible with respect to their values in the slipping modes. Therefore, the region of phase space including the states in which the asperities are stationary (sticking region) is a subset of the space $XYZ$. It is the region bounded by the planes $X$ $=$ 0, $Y$ $=$ 0, ${\mathrm{\Pi }}_{\mathrm{1}}$ and ${\mathrm{\Pi }}_{\mathrm{2}}$: a tetrahedron $T$ (Fig. 2).

A seismic event takes place when the orbit of the system reaches one of the faces $ACD$ or $BCD$ of $T$, belonging, respectively, to the planes ${\mathrm{\Pi }}_{\mathrm{1}}$ and ${\mathrm{\Pi }}_{\mathrm{2}}$. In these cases, the system passes from mode 00 to mode 10 or 01, respectively. If the orbit reaches the edge $CD$, the system passes instead to mode 11. For later use, we introduce a point $P$ with coordinates $$X_{\text{P}}=\frac{\mathit{\alpha }+\mathit{\alpha }\mathit{\beta }+\mathrm{1}}{\mathrm{1}+\mathrm{2}\mathit{\alpha }},Y_{\text{P}}=\frac{\mathit{\alpha }+\mathit{\alpha }\mathit{\beta }+\mathit{\beta }}{\mathrm{1}+\mathrm{2}\mathit{\alpha }},Z_{\text{P}}=-\frac{\mathrm{1}-\mathit{\beta }}{\mathrm{1}+\mathrm{2}\mathit{\alpha }}.$$ It belongs to the edge $CD$ and corresponds to the case of elastic coupling; in fact $Z_{\text{P}}$ $=$ $Y_{\text{P}}$ $-$ $X_{\text{P}}$.

Equations of motion and solutions

The equations of motions of the four dynamic modes and the corresponding solutions are given below. Viscoelastic relaxation is negligible during the slipping modes; therefore, the equations for $X$ and $Y$ are the same as in the case of elastic coupling, while $Z$ changes according to the equation $\ddot{Z}$ $=$ $\ddot{Y}$ $-$ $\ddot{X}$.

The sticking region $\mathbf{T}$ of the system is a tetrahedron $ABCD$ in the phase space ($\mathit{\alpha }$ $=$ 1, $\mathit{\beta }$ $=$ 1). The point $P$ is indicated.

[Figure omitted. See PDF]

Stationary asperities (mode 00)

The variables $X$ and $Y$ increase steadily due to tectonic motions, while $Z$ is governed by the Maxwell constitutive equation. The equations of motion are $$\ddot{X}=\mathrm{0},\ddot{Y}=\mathrm{0},\ddot{Z}=\frac{Z}{{\mathrm{\Theta }}^{\mathrm{2}}},$$ where a dot indicates differentiation with respect to $T$. The fault can enter mode 00 from mode 10 or from mode 01. With initial conditions $$\begin{array}{}& & X\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{X},Y\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{Y},Z\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{Z},& & \dot{X}\left(\mathrm{0}\right)=V,\dot{Y}\left(\mathrm{0}\right)=V,\dot{Z}\left(\mathrm{0}\right)=-\frac{\stackrel{\mathrm{‾}}{Z}}{\mathrm{\Theta }},\end{array}$$ the solution is $$X\left(T\right)=\stackrel{\mathrm{‾}}{X}+VT,Y\left(T\right)=\stackrel{\mathrm{‾}}{Y}+VT,Z\left(T\right)=\stackrel{\mathrm{‾}}{Z}{e}^{-T/\mathrm{\Theta }},$$ with $T$ $\ge$ 0. The initial point belongs necessarily to $T$ and Eq. () are the parametric equations of a curve lying on the plane $$X-Y+\stackrel{\mathrm{‾}}{Y}-\stackrel{\mathrm{‾}}{X}=\mathrm{0},$$ which is parallel to the $Z$ axis.

Failure of asperity 1 (mode 10)

The equations of motion are $$\begin{array}{}& & \ddot{X}+X-\mathit{\alpha }Z-\mathit{ϵ}=\mathrm{0},& & \ddot{Y}=\mathrm{0},& & \ddot{Z}-X+\mathit{\alpha }Z+\mathit{ϵ}=0.\end{array}$$ The fault can enter mode 10 from mode 11 or from mode 00.

• a.

  Case 11 $\to$ 10: with initial conditions$$\begin{array}{}& & X\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{X},Y\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{Y},Z\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{Z},& & \dot{X}\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{V},\dot{Y}\left(\mathrm{0}\right)=\mathrm{0},\dot{Z}\left(\mathrm{0}\right)=-\stackrel{\mathrm{‾}}{V},\end{array}$$

  the solution is$$\begin{array}{}& & X\left(T\right)=\stackrel{\mathrm{‾}}{X}-\frac{{\stackrel{\mathrm{‾}}{U}}_{\mathrm{1}}}{\mathrm{2}}(\mathrm{1}-\cos \mathit{\omega }T)+\frac{\stackrel{\mathrm{‾}}{V}}{\mathit{\omega }}\sin \mathit{\omega }T,& & Y\left(T\right)=\stackrel{\mathrm{‾}}{Y},& & Z\left(T\right)=\stackrel{\mathrm{‾}}{Z}+\stackrel{\mathrm{‾}}{X}-X\left(T\right),\end{array}$$

  where$${\stackrel{\mathrm{‾}}{U}}_{\mathrm{1}}=\mathrm{2}\frac{\stackrel{\mathrm{‾}}{X}-\mathit{\alpha }\stackrel{\mathrm{‾}}{Z}-\mathit{ϵ}}{{\mathit{\omega }}^{\mathrm{2}}}.$$

  The slip duration, calculated from the condition $\dot{X}\left(T\right)$ $=$ 0, is$$T_{10}=\frac{\mathrm{1}}{\mathit{\omega }}\left(\mathit{\pi },+,\arctan ,\frac{\mathrm{2}\stackrel{\mathrm{‾}}{V}}{\mathit{\omega }{\stackrel{\mathrm{‾}}{U}}_{\mathrm{1}}}\right),$$

  and the final slip amplitude is$$U_{10}=\frac{{\stackrel{\mathrm{‾}}{U}}_{\mathrm{1}}}{\mathrm{2}}+\sqrt{{\left(\frac{{\stackrel{\mathrm{‾}}{U}}_{\mathrm{1}}}{\mathrm{2}}\right)}^{\mathrm{2}}+{\left(\frac{\stackrel{\mathrm{‾}}{V}}{\mathit{\omega }}\right)}^{\mathrm{2}}}.$$

• b.

  Case 00 $\to$ 10: in this case the initial point belongs to the face $ACD$ so that$$\stackrel{\mathrm{‾}}{X}-\mathit{\alpha }\stackrel{\mathrm{‾}}{Z}=\mathrm{1},\stackrel{\mathrm{‾}}{V}=\mathrm{0},$$

  and from Eq. ()$${\stackrel{\mathrm{‾}}{U}}_{\mathrm{1}}=U.$$

  The solution reduces to$$\begin{array}{}& & X\left(T\right)=\stackrel{\mathrm{‾}}{X}-\frac{U}{\mathrm{2}}(\mathrm{1}-\cos \mathit{\omega }T),& & Y\left(T\right)=\stackrel{\mathrm{‾}}{Y},& & Z\left(T\right)=\stackrel{\mathrm{‾}}{Z}+\frac{U}{\mathrm{2}}(\mathrm{1}-\cos \mathit{\omega }T).\end{array}$$

  If the orbit does not reach the face $BCD$ during the mode, one has$$T_{10}=\frac{\mathit{\pi }}{\mathit{\omega }},U_{10}=U.$$

  If the orbit reaches $BCD$ before time $\mathit{\pi }/\mathit{\omega }$ has elapsed, the system passes to mode 11. In this case,$$T_{10}=\frac{\mathit{\pi }}{\mathit{\omega }}-\frac{\mathrm{1}}{\mathit{\omega }}\arccos \left(\mathrm{2},\frac{U_{10}}{U},-,\mathrm{1}\right),$$

  where$$U_{10}=\frac{\mathit{\beta }-\stackrel{\mathrm{‾}}{Y}-\mathit{\alpha }\stackrel{\mathrm{‾}}{Z}}{\mathit{\alpha }}.$$

Failure of asperity 2 (mode 01)

The equations of motion are $$\begin{array}{}& & \ddot{X}=\mathrm{0},& & \ddot{Y}+Y+\mathit{\alpha }Z-\mathit{\beta }\mathit{ϵ}=\mathrm{0},& & \ddot{Z}+Y+\mathit{\alpha }Z-\mathit{\beta }\mathit{ϵ}=0.\end{array}$$ The fault can enter mode 01 from mode 11 or from mode 00.

• a.

  Case 11 $\to$ 01: with initial conditions $$\begin{array}{}& & X\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{X},Y\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{Y},Z\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{Z},& & \dot{X}\left(\mathrm{0}\right)=\mathrm{0},\dot{Y}\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{V},\dot{Z}\left(\mathrm{0}\right)=\stackrel{\mathrm{‾}}{V},\end{array}$$

  the solution is $$\begin{array}{}& & X\left(T\right)=\stackrel{\mathrm{‾}}{X},& & Y\left(T\right)=\stackrel{\mathrm{‾}}{Y}-\frac{{\stackrel{\mathrm{‾}}{U}}_{\mathrm{2}}}{\mathrm{2}}(\mathrm{1}-\cos \mathit{\omega }T)+\frac{\stackrel{\mathrm{‾}}{V}}{\mathit{\omega }}\sin \mathit{\omega }T,& & Z\left(T\right)=\stackrel{\mathrm{‾}}{Z}-\stackrel{\mathrm{‾}}{Y}+Y\left(T\right),\end{array}$$

  where $${\stackrel{\mathrm{‾}}{U}}_{\mathrm{2}}=\mathrm{2}\frac{\stackrel{\mathrm{‾}}{Y}+\mathit{\alpha }\stackrel{\mathrm{‾}}{Z}-\mathit{\beta }\mathit{ϵ}}{{\mathit{\omega }}^{\mathrm{2}}}.$$

  The slip duration, calculated from the condition $\dot{Y}\left(T\right)$ $=$ 0, is $$T_{01}=\frac{\mathrm{1}}{\mathit{\omega }}\left(\mathit{\pi },+,\arctan ,\frac{\mathrm{2}\stackrel{\mathrm{‾}}{V}}{\mathit{\omega }{\stackrel{\mathrm{‾}}{U}}_{\mathrm{2}}}\right),$$

  and the final slip amplitude is $$U_{01}=\frac{{\stackrel{\mathrm{‾}}{U}}_{\mathrm{2}}}{\mathrm{2}}+\sqrt{{\left(\frac{{\stackrel{\mathrm{‾}}{U}}_{\mathrm{2}}}{\mathrm{2}}\right)}^{\mathrm{2}}+{\left(\frac{\stackrel{\mathrm{‾}}{V}}{\mathit{\omega }}\right)}^{\mathrm{2}}}.$$

• b.

  Case 00 $\to$ 01: in this case the initial point belongs to the face $BCD$ so that $$\stackrel{\mathrm{‾}}{Y}+\mathit{\alpha }\stackrel{\mathrm{‾}}{Z}=\mathit{\beta },\stackrel{\mathrm{‾}}{V}=\mathrm{0},$$

  and from Eq. () $${\stackrel{\mathrm{‾}}{U}}_{\mathrm{2}}=\mathit{\beta }U.$$

  The solution reduces to $$\begin{array}{}& & X\left(T\right)=\stackrel{\mathrm{‾}}{X},& & Y\left(T\right)=\stackrel{\mathrm{‾}}{Y}-\frac{\mathit{\beta }U}{\mathrm{2}}(\mathrm{1}-\cos \mathit{\omega }T),& & Z\left(T\right)=\stackrel{\mathrm{‾}}{Z}-\frac{\mathit{\beta }U}{\mathrm{2}}(\mathrm{1}-\cos \mathit{\omega }T).\end{array}$$

  If the orbit does not reach the face $ACD$ during the mode, one has $$T_{01}=\frac{\mathit{\pi }}{\mathit{\omega }},U_{01}=\mathit{\beta }U.$$

  If the orbit reaches $ACD$ before time $\mathit{\pi }/\mathit{\omega }$ has elapsed, the system passes to mode 11. In this case, $$T_{01}=\frac{\mathit{\pi }}{\mathit{\omega }}-\frac{\mathrm{1}}{\mathit{\omega }}\arccos \left(\mathrm{2},\frac{U_{01}}{\mathit{\beta }U},-,\mathrm{1}\right),$$

  where $$U_{01}=\frac{\mathrm{1}-\stackrel{\mathrm{‾}}{X}+\mathit{\alpha }\stackrel{\mathrm{‾}}{Z}}{\mathit{\alpha }}.$$

Simultaneous asperity failure (mode 11)

The equations of motion are $$\begin{array}{}& & \ddot{X}+X-\mathit{\alpha }Z-\mathit{ϵ}=\mathrm{0},& & \ddot{Y}+Y+\mathit{\alpha }Z-\mathit{\beta }\mathit{ϵ}=\mathrm{0},& & \ddot{Z}-X+Y+\mathrm{2}\mathit{\alpha }Z+(\mathrm{1}-\mathit{\beta })\mathit{ϵ}=\mathrm{0},\end{array}$$ and the solution is $$\begin{array}{}& & X\left(T\right)=A\sin T+B\cos T+C\sin \mathrm{\Omega }T+D\cos \mathrm{\Omega }T+E_{\mathrm{1}},& & Y\left(T\right)=A\sin T+B\cos T-C\sin \mathrm{\Omega }T-D\cos \mathrm{\Omega }T+E_{\mathrm{2}},& & Z\left(T\right)=-\mathrm{2}C\sin \mathrm{\Omega }T-\mathrm{2}D\cos \mathrm{\Omega }T+E_{\mathrm{3}},\end{array}$$ where the constants $A$, $B$, $C$, $D$, $E_{\mathrm{1}}$, $E_{\mathrm{2}}$, and $E_{\mathrm{3}}$ depend on initial conditions. The fault can enter mode 11 from mode 10, 01 or 00.

• a.

  Case 10 $\to$ 11: the initial conditions are $$\begin{array}{}& & X=\stackrel{\mathrm{‾}}{X},Y=\stackrel{\mathrm{‾}}{Y},Z=\stackrel{\mathrm{‾}}{Z},& & \dot{X}=\stackrel{\mathrm{‾}}{V},\dot{Y}=\mathrm{0},\dot{Z}=-\stackrel{\mathrm{‾}}{V},\end{array}$$

  and the constants are $$\begin{array}{}& & B=\frac{\mathrm{1}}{\mathrm{2}}\left[\stackrel{\mathrm{‾}}{X},+,\stackrel{\mathrm{‾}}{Y},-,\mathit{ϵ},\left(X_{\text{P}},+,Y_{\text{P}}\right)\right],& & D=\frac{\mathrm{1}}{\mathrm{2}}\left(\mathit{ϵ},Z_{\text{P}},+,\frac{\stackrel{\mathrm{‾}}{X}-\stackrel{\mathrm{‾}}{Y}-\mathrm{2}\mathit{\alpha }\stackrel{\mathrm{‾}}{Z}}{{\mathrm{\Omega }}^{\mathrm{2}}}\right),& & E_{\mathrm{1}}=\mathit{ϵ}{X}_{\text{P}}+\mathit{\alpha }\frac{\stackrel{\mathrm{‾}}{X}-\stackrel{\mathrm{‾}}{Y}+\stackrel{\mathrm{‾}}{Z}}{{\mathrm{\Omega }}^{\mathrm{2}}},& & E_{\mathrm{2}}=\mathit{ϵ}{Y}_{\text{P}}-\mathit{\alpha }\frac{\stackrel{\mathrm{‾}}{X}-\stackrel{\mathrm{‾}}{Y}+\stackrel{\mathrm{‾}}{Z}}{{\mathrm{\Omega }}^{\mathrm{2}}},& & E_{\mathrm{3}}=\mathit{ϵ}{Z}_{\text{P}}+\frac{\stackrel{\mathrm{‾}}{X}-\stackrel{\mathrm{‾}}{Y}+\stackrel{\mathrm{‾}}{Z}}{{\mathrm{\Omega }}^{\mathrm{2}}},& & A=\frac{\stackrel{\mathrm{‾}}{V}}{\mathrm{2}},C=\frac{\stackrel{\mathrm{‾}}{V}}{\mathrm{2}\mathrm{\Omega }}.\end{array}$$

• b.

  Case 01 $\to$ 11: the initial conditions are $$\begin{array}{}& & X=\stackrel{\mathrm{‾}}{X},Y=\stackrel{\mathrm{‾}}{Y},Z=\stackrel{\mathrm{‾}}{Z},& & \dot{X}=\mathrm{0},\dot{Y}=\stackrel{\mathrm{‾}}{V},\dot{Z}=\stackrel{\mathrm{‾}}{V}.\end{array}$$

  The constants $B$, $D$, $E_{\mathrm{1}}$, $E_{\mathrm{2}}$, and $E_{\mathrm{3}}$ are given by Eqs. ()–(), while $$A=\frac{\stackrel{\mathrm{‾}}{V}}{\mathrm{2}},C=-\frac{\stackrel{\mathrm{‾}}{V}}{\mathrm{2}\mathrm{\Omega }}.$$

• c.

  Case 00 $\to$ 11: the initial conditions are $$\begin{array}{}& & X=\stackrel{\mathrm{‾}}{X},Y=\stackrel{\mathrm{‾}}{Y},Z=\stackrel{\mathrm{‾}}{Z},& & \dot{X}=\mathrm{0},\dot{Y}=\mathrm{0},\dot{Z}=0.\end{array}$$

  The constants $B$, $D$, $E_{\mathrm{1}}$, $E_{\mathrm{2}}$, and $E_{\mathrm{3}}$ are given by Eqs. ()–(), while $$A=\mathrm{0},C=0.$$

The sequence of slipping modes

In general, a seismic event will involve $n$ slipping modes of the fault. The sequence of slipping modes determines not only the source function and the seismic moment of the earthquake, but also the position of its focus. We wish to find the relationship between the state of the fault before the earthquake and the sequence of slipping modes.

The surface $\mathrm{\Sigma }$ divides the sticking region $\mathbf{T}$ in two subsets ${\mathbf{T}}_{\mathrm{1}}$ (below) and ${\mathbf{T}}_{\mathrm{2}}$ (above), which determine the first slipping mode of the seismic event ($\mathit{\alpha }$ $=$ 1, $\mathit{\beta }$ $=$ 1).

[Figure omitted. See PDF]

During the interseismic intervals, the fault is subject to continuous tectonic loading due to the motion of adjacent plates and to the effect of viscoelastic relaxation of the stress accumulated by previous seismic activity. Given any state $P_{\mathrm{0}}$ $=$ ($X_{\mathrm{0}}$, $Y_{\mathrm{0}}$, $Z_{\mathrm{0}}$) $\in$ $\mathbf{T}$, its orbit will lead to the failure of asperity 1 or asperity 2 or to the simultaneous failure of both asperities. In fact, all the orbits (Eq. ) in mode 00 reach one of the faces $ACD$ or $BCD$ or their common edge $CD$. We wish to determine the subset ${\mathbf{T}}_{\mathrm{1}}$ of the sticking region $\mathbf{T}$ such that, if $P_{\mathrm{0}}$ $\in$ ${\mathbf{T}}_{\mathrm{1}}$, the orbit reaches $ACD$ and the subset ${\mathbf{T}}_{\mathrm{2}}$ such that, if $P_{\mathrm{0}}$ $\in$ ${\mathbf{T}}_{\mathrm{2}}$, the orbit reaches $BCD$.

Any curve (Eq. ), if prolonged outside $\mathbf{T}$, intersects both ${\mathrm{\Pi }}_{\mathrm{1}}$ and ${\mathrm{\Pi }}_{\mathrm{2}}$. Let $P_{\mathrm{1}}$ $=$ ($X_{\mathrm{1}}$, $Y_{\mathrm{1}}$, $Z_{\mathrm{1}}$) and $P_{\mathrm{2}}$ $=$ ($X_{\mathrm{2}}$, $Y_{\mathrm{2}}$, $Z_{\mathrm{2}}$) be the intersection points with the two planes, respectively, and let $T_{\mathrm{1}}$ and $T_{\mathrm{2}}$ be the corresponding instants in time. Accordingly, $X_{\mathrm{1}}$ and $Z_{\mathrm{1}}$ must satisfy Eq. () or, thanks to Eq. (), $$X_{\mathrm{0}}+V{T}_{\mathrm{1}}=\mathrm{1}+\mathit{\alpha }{Z}_{\mathrm{0}}{e}^{-T_{\mathrm{1}}/\mathrm{\Theta }},$$ whence $$T_{\mathrm{1}}=\mathrm{\Theta }W\left({\mathit{\gamma }}_{\mathrm{1}}\right)+\frac{\mathrm{1}-X_{\mathrm{0}}}{V},$$ where $W$ is the Lambert function with argument $${\mathit{\gamma }}_{\mathrm{1}}=\frac{\mathit{\alpha }{Z}_{\mathrm{0}}}{V\mathrm{\Theta }}{e}^{-\frac{\mathrm{1}-X_{\mathrm{0}}}{V\mathrm{\Theta }}}.$$ Analogously, $Y_{\mathrm{2}}$ and $Z_{\mathrm{2}}$ must satisfy Eq. () or, thanks to Eq. (), $$Y_{\mathrm{0}}+V{T}_{\mathrm{2}}=\mathit{\beta }-\mathit{\alpha }{Z}_{\mathrm{0}}{e}^{-T_{\mathrm{2}}/\mathrm{\Theta }},$$ whence $$T_{\mathrm{2}}=\mathrm{\Theta }W\left({\mathit{\gamma }}_{\mathrm{2}}\right)+\frac{\mathit{\beta }-Y_{\mathrm{0}}}{V}$$ with $${\mathit{\gamma }}_{\mathrm{2}}=-\frac{\mathit{\alpha }{Z}_{\mathrm{0}}}{V\mathrm{\Theta }}{e}^{-\frac{\mathit{\beta }-Y_{\mathrm{0}}}{V\mathrm{\Theta }}}.$$ We consider the difference $$\mathrm{\Delta }T=T_{\mathrm{1}}-T_{\mathrm{2}}$$ and define a surface $\mathrm{\Sigma }$ with the equation $$\mathrm{\Delta }T(X,Y,Z)=\mathrm{0}$$ or, thanks to Eqs. () and (), $$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 }=0.$$ This is a transcendental surface that divides $\mathbf{T}$ in two connected, open subsets ${\mathbf{T}}_{\mathrm{1}}$ and ${\mathbf{T}}_{\mathrm{2}}$ with the required properties (Fig. 3). If $\mathit{\beta }$ $=$ 1, the surface $\mathrm{\Sigma }$ divides $\mathbf{T}$ in two halves; if $\mathit{\beta }$ $<$ 1, ${\mathbf{T}}_{\mathrm{1}}$ has a smaller volume than ${\mathbf{T}}_{\mathrm{2}}$. The edge $CD$ of $\mathbf{T}$ belongs to $\mathrm{\Sigma }$. By definition, no orbit can cross $\mathrm{\Sigma }$; therefore, if $P_{\mathrm{0}}$ $\in$ $\mathrm{\Sigma }$, its orbit remains on $\mathrm{\Sigma }$ and reaches the edge $CD$.

The faces $ACD$ and $BCD$ of $\mathbf{T}$ and their subsets, which determine the sequence of slipping modes and the moment rate of the seismic event ($\mathit{\alpha }$ $=$ 1, $\mathit{\beta }$ $=$ 1, $\mathit{ϵ}$ $=$ 0.7).

[Figure omitted. See PDF]

After an orbit reaches one of the faces $ACD$ or $BCD$ at a point $P_k$, the sequence of modes in the earthquake will be different according to which subset of the face $P_k$ belongs to. This is illustrated in Fig. 4. Let us consider the face $ACD$. If $P_k$ belongs to the triangle $Q_{\mathrm{1}}$, the earthquake will be a one-mode event 10. If $P_k$ belongs to the segment $s_{\mathrm{1}}$, the earthquake will be a two-mode event 10–01. If $P_k$ belongs to the trapezoid $R_{\mathrm{1}}$, the earthquake will be a three-mode event 10–11–10 or 10–11–01. The remaining part of the face would lead to overshooting. Analogous considerations can be made for subsets $Q_{\mathrm{2}}$, $s_{\mathrm{2}}$ and $R_{\mathrm{2}}$ of the face $BCD$.

Therefore, events involving the simultaneous failure of asperities can take place only from particular subsets of states of the system. In general, a three-mode event can result from four different sequences of modes: 10–11–10, 10–11–01, 01–11–10, and 01–11–01. In particular cases, sequences like 10–01–10 or 01–10–01 are also possible.

The reasons for the different sequences of modes involved in the earthquake are clear if we consider the forces acting on the asperities in the different states. If we consider the face $ACD$, we have $F_{\mathrm{1}}$ $=$ $-$1 everywhere, while $F_{\mathrm{2}}$ is equal to $-\mathit{\beta }$ on $CD$ and decreases in magnitude with the distance from this edge. Hence, the onset mode 10 of the sequence can trigger mode 11 only if $\left|F_{\mathrm{2}}\right|$ is large enough and this occurs if $P_k$ $\in$ $R_{\mathrm{1}}$. When $F_{\mathrm{2}}$ $=$ $-(\mathit{\beta }$ $-$ $\mathit{\alpha }U)$ we have the limit case of two consecutive modes 10–01. For smaller values of $\left|F_{\mathrm{2}}\right|$, no triggering occurs and the earthquake is a one-mode event 10. The same considerations can be made for the face $BCD$.

This analysis enlightens the relationship between the state of the fault before an earthquake and the sequence of modes in the seismic event. It also suggests that the knowledge of the source function of an earthquake may allow us to constrain the orbit of the system in the phase space.

Seismic moment rates

The number and the sequence of slipping modes involved in a seismic event determine the moment rate of the earthquake. Let $P_i$ be the singular points of the orbit, i.e. the points where the system passes from one mode to another. If the seismic event begins at $P_k$, the representative point of the system when it enters the $i$th slipping mode is $P_{k+i-\mathrm{1}}$ and the corresponding instant in time is $T_{k+i-\mathrm{1}}$ ($i$ $=$ 1, 2, … $n$). The duration of the $i$th mode is $$\mathrm{\Delta }{T}_i=T_{k+i}-T_{k+i-\mathrm{1}},$$ and the seismic event terminates at time $T_{k+n}$. In the $i$th mode, the slip functions of asperities 1 and 2 are, respectively, $$\begin{array}{}& & \mathrm{\Delta }{X}_i\left(T\right)=X_{k+i-\mathrm{1}}-X\left(T,-,T_{k+i-\mathrm{1}}\right),& & \mathrm{\Delta }{Y}_i\left(T\right)=Y_{k+i-\mathrm{1}}-Y\left(T,-,T_{k+i-\mathrm{1}}\right),\end{array}$$ where the appropriate expressions of $X\left(T\right)$ and $Y\left(T\right)$ must be used. The moment rate of an $n$-mode seismic event can be calculated as $$\begin{array}{cc}\dot{M}\left(T\right)=& \frac{M_{\mathrm{1}}}{U}\sum _{i=\mathrm{1}}^n\left(\mathrm{\Delta },{\dot{X}}_i,+,\mathrm{\Delta },{\dot{Y}}_i\right)\\ & & \left[H,\left(T,-,T_{k+i-\mathrm{1}}\right),-,H,\left(T,-,T_{k+i}\right)\right],\end{array}$$ where $M_{\mathrm{1}}$ is the seismic moment due to the slip of asperity 1 by an amount $U$ and $H\left(T\right)$ is the Heaviside function. The final slip amplitudes of asperities 1 and 2 are, respectively, $$U_{\mathrm{1}}=\sum _{i=\mathrm{1}}^n\mathrm{\Delta }{X}_i\left(T_{k+i}\right),U_{\mathrm{2}}=\sum _{i=\mathrm{1}}^n\mathrm{\Delta }{Y}_i\left(T_{k+i}\right),$$ and the final seismic moment is $$M_{\mathrm{0}}=M_{\mathrm{1}}\frac{U_{\mathrm{1}}+U_{\mathrm{2}}}{U}.$$ The moment rate depends on the state of the fault at the beginning of the seismic event, i.e. on the coordinates $X_k$, $Y_k$ and $Z_k$. This state is a priori unknown, but the knowledge of the source function of the earthquake allows us to set constraints on it. As shown in Sect. 4, if the first mode is 10 or 01, $P_k$ must belong to the face $ACD$ or $BCD$ of $\mathbf{T}$. In addition, if the event has a single mode, $P_k$ belongs to the subset $Q_{\mathrm{1}}$ or $Q_{\mathrm{2}}$; if the event has two modes, $P_k$ belongs to the segment $s_{\mathrm{1}}$ or $s_{\mathrm{2}}$; if the event has three modes, $P_k$ belongs to the subset $R_{\mathrm{1}}$ or $R_{\mathrm{2}}$.

This allows us to constrain the evolution of the system to a certain subset of the phase space and, when the next earthquake occurs, the knowledge of its moment rate will allow us to further constrain this subset. Hence, if we knew the source functions of a sufficiently large number of consecutive earthquakes, we could constrain more and more the orbit of the system and its evolution could be predicted with a smaller uncertainty.

Application to the 1964 Alaska earthquake

The 1964 Alaska earthquake was one of the largest earthquakes in the last century, with a seismic moment $M_{\mathrm{0}}$ $=$ 3 $\times$ 10${}^{22}$ Nm (Christensen and Beck, 1994; Holdahl and Sauber, 1994; Johnson et al., 1996; Ichinose et al., 2007). Seismological, geodetic and tsunami data indicate that the earthquake was the result of the slipping of two asperities, the Prince William Sound and the Kodiak Island asperity, which we call asperities 1 and 2, respectively. The earthquake started with the failure of asperity 1 followed by that of asperity 2. On the basis of coseismic surface deformation, Santini et al. (2003) suggested average slips of $u_{\mathrm{1}}$ $=$ 24 m for asperity 1 and $u_{\mathrm{2}}$ $=$ 18 m for asperity 2.

For the Alaska earthquake there is clear evidence that post-seismic deformation occurred in the decades following the event (Zweck et al., 2002; Suito and Freymueller, 2009). Part of the deformation has been ascribed to aseismic slip of the fault and part to viscoelastic relaxation. The latter shows a characteristic time $\mathit{\tau }$ $\approx$ 30 a. The relative plate velocity is $v$ $=$ 5.7 cm a${}^{-\mathrm{1}}$ (DeMets and Dixon, 1999; Cohen and Freymueller, 2004). In fact, the velocity of the Pacific Plate relative to the North American Plate at the Alaska/Aleutian Trench increases gradually from the northeast to the southwest. However, the difference between the area of Prince William Sound and the area of Kodiak Island is small, in the order of few millimetres per year, and can be reasonably neglected.

According to the present model, the seismic event was a sequence of modes 10–01 starting from mode 00. Since the first mode was 10, the orbit of the system in mode 00 was in the subset ${\mathbf{T}}_{\mathrm{1}}$ of the sticking region. Let $P_{\mathrm{1}}$ be the representative point of the system at the beginning of the seismic event. Since mode 10 was followed by mode 01, $P_{\mathrm{1}}$ belongs to segment $s_{\mathrm{1}}$ (Fig. 4). We may express the coordinates of $P_{\mathrm{1}}$ as $$X_{\mathrm{1}}=\mathit{\alpha }{Z}_{\mathrm{1}}+\mathrm{1},Y_{\mathrm{1}}=\mathit{\beta }-\mathit{\alpha }U-\mathit{\alpha }{Z}_{\mathrm{1}},Z_{\mathrm{1}}$$ with $$Z_{\text{a}}\le Z_{\mathrm{1}}\le Z_{\text{b}},$$ where $$Z_{\text{a}}=-\frac{\mathrm{1}-U}{\mathit{\alpha }},Z_{\text{b}}=\frac{\mathit{\beta }-(\mathit{\alpha }+\mathit{\beta })U}{\mathit{\alpha }}.$$ The orbit of the system is one curve of the bundle with parametric equations (Eq. ) passing through $s_{\mathrm{1}}$. At the end of mode 10, the system is at $P_{\mathrm{2}}$ with coordinates $$X_{\mathrm{2}}=\mathit{\alpha }{Z}_{\mathrm{1}}+\mathrm{1}-U,Y_{\mathrm{2}}=\mathit{\beta }-\mathit{\alpha }U-\mathit{\alpha }{Z}_{\mathrm{1}},Z_{\mathrm{2}}=Z_{\mathrm{1}}+U.$$ As $Z_{\mathrm{1}}$ varies in the interval Eq. (), there is an infinite number of points $P_{\mathrm{2}}$ forming another segment $r_{\mathrm{1}}$ belonging to the face $BCD$ and parallel to the edge $CD$. At the end of the event, the system is at $P_{\mathrm{3}}$, with coordinates $$\begin{array}{cc}{X}_{\mathrm{3}}=& \mathit{\alpha }{Z}_{\mathrm{1}}+\mathrm{1}-U,Y_{\mathrm{3}}=\mathit{\beta }-(\mathit{\alpha }+\mathit{\beta })U-\mathit{\alpha }{Z}_{\mathrm{1}},\\ & & Z_{\mathrm{3}}=Z_{\mathrm{1}}+(\mathrm{1}-\mathit{\beta })U.\end{array}$$ As $Z_{\mathrm{1}}$ varies in the interval Eq. (), there is an infinite number of points $P_{\mathrm{3}}$ forming another segment $q_{\mathrm{1}}$. This segment is also parallel to the edge $CD$. However, it intersects the surface $\mathrm{\Sigma }$ for $Z_{\mathrm{1}}$ $=$ $Z_{\text{c}}$, with $Z_{\text{a}}$ $<$ $Z_{\text{c}}$ $<$ $Z_{\text{b}}$.

From Eq. (), it is easy to calculate the forces on the asperities at points $P_{\mathrm{1}}$, $P_{\mathrm{2}}$ and $P_{\mathrm{3}}$. These forces are independent of the positions of the points on the respective segments $s_{\mathrm{1}}$, $r_{\mathrm{1}}$ and $q_{\mathrm{1}}$: $$\begin{array}{}& & F_{\mathrm{1}}=-\mathrm{1},F_{\mathrm{2}}=-(\mathit{\beta }-\mathit{\alpha }U)\text{on}{s}_{\mathrm{1}},& & F_{\mathrm{1}}=-(\mathrm{2}\mathit{ϵ}-\mathrm{1}),F_{\mathrm{2}}=-\mathit{\beta }\text{on}{r}_{\mathrm{1}},& & F_{\mathrm{1}}=-(\mathrm{2}\mathit{ϵ}-\mathrm{1}+\mathit{\alpha }\mathit{\beta }U),F_{\mathrm{2}}=-(\mathrm{2}\mathit{ϵ}-\mathrm{1})\mathit{\beta }\text{on}{q}_{\mathrm{1}}.\end{array}$$ For an application of the model to the Alaska earthquake, we take $\mathit{\alpha }$ $=$ 0.1, $\mathit{\beta }$ $=$ 0.75, and $\mathit{ϵ}$ $=$ 0.7 (Dragoni and Santini, 2012). It follows $U$ $\simeq$ 0.545 and $V\mathrm{\Theta }$ $\simeq$ 0.039. With these values, Eq. () yields $Z_{\text{a}}$ $\simeq$ $-$4.55 and $Z_{\text{b}}$ $\simeq$ 2.86, while $Z_{\text{c}}$ $\simeq$ 0.41.

Then, according to Eqs. ()–(), the forces immediately before the 1964 earthquake are $F_{\mathrm{1}}\left(T_{\mathrm{1}}\right)$ $=$ $-$1 and $F_{\mathrm{2}}\left(T_{\mathrm{1}}\right)$ $=$ $-$0.70, showing that the magnitude of stress on asperity 2 is 70 % of that on asperity 1. The failure of asperity 1 reduces the stress on asperity 1 and transfers stress to asperity 2 up to static friction, so that $F_{\mathrm{1}}\left(T_{\mathrm{2}}\right)$ $=$ $-$0.40 and $F_{\mathrm{2}}\left(T_{\mathrm{2}}\right)$ $=$ $-$0.75. Finally, the failure of asperity 2 reduces the stress on asperity 2 and transfers stress back to asperity 1, so that at the end of the event it results in $F_{\mathrm{1}}\left(T_{\mathrm{3}}\right)$ $=$ $-$0.44 and $F_{\mathrm{2}}\left(T_{\mathrm{3}}\right)$ $=$ $-$0.30, indicating a more homogeneous stress distribution. Then the system evolves in mode 00, where both stresses increase in magnitude, but at different rates.

Post-seismic evolution

On the basis of a purely elastic model, Dragoni and Santini (2012) predicted that the next large earthquake involving the 1964 fault would take place about 350 years later and would be due to the failure of asperity 2 alone. If we introduce viscoelastic relaxation, a wider range of possibilities appears. Since the segment $q_{\mathrm{1}}$ intersects $\mathrm{\Sigma }$, the point $P_{\mathrm{3}}$ can belong to ${\mathbf{T}}_{\mathrm{1}}$, ${\mathbf{T}}_{\mathrm{2}}$ or $\mathrm{\Sigma }$. 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 both asperities.

(a) Duration $T^{\prime }$ of the interseismic interval following an event with mode sequence 10–01, (b) forces $F_{\mathrm{1}}$ and $F_{\mathrm{2}}$ on the asperities at the beginning of the subsequent event, and (c) seismic moment $M_{\mathrm{0}}$ of the subsequent event, as functions of the variable $Z_{\mathrm{1}}$ characterizing the initial state of the system. The values of parameters are appropriate to the 1964 Alaska earthquake ($\mathit{\alpha }$ $=$ 0.1, $\mathit{\beta }$ $=$ 0.75, $\mathit{ϵ}$ $=$ 0.7, $V\mathrm{\Theta }$ $=$ 0.039).

[Figure omitted. See PDF]

Examples of possible moment rates $\dot{M}\left(T\right)$ for the event following the 10–01 event: (a) $-$4.55 $\le$ $Z_{\mathrm{1}}$ $<$ 0.20, (b) $Z_{\mathrm{1}}$ $=$ 0.30, (c) $Z_{\mathrm{1}}$ $=$ 0.41, (d) $Z_{\mathrm{1}}$ $=$ 0.60, and (e) 0.70 $\le$ $Z_{\mathrm{1}}$ $\le$ 2.86, where $Z_{\mathrm{1}}$ characterizes the initial state of the system. The values of parameters are appropriate to the 1964 Alaska earthquake ($\mathit{\alpha }$ $=$ 0.1, $\mathit{\beta }$ $=$ 0.75, $\mathit{ϵ}$ $=$ 0.7, $V\mathrm{\Theta }$ $=$ 0.039).

[Figure omitted. See PDF]

According to the present model, the duration of the interseismic interval between 1964 and the next earthquake is $$\frac{T^{\prime }}{\mathrm{\Theta }}=\left\{\begin{array}{ll}W\left({\mathit{\gamma }}_{\mathrm{1}}^{\prime }\right)+\frac{\mathrm{1}-X_{\mathrm{3}}}{V\mathrm{\Theta }},& P_{\mathrm{3}}\in {\mathbf{T}}_{\mathrm{1}}\\ W\left({\mathit{\gamma }}_{\mathrm{2}}^{\prime }\right)+\frac{\mathit{\beta }-Y_{\mathrm{3}}}{V\mathrm{\Theta }},& P_{\mathrm{3}}\in {\mathbf{T}}_{\mathrm{2}}\end{array},,\right.$$ where $${\mathit{\gamma }}_{\mathrm{1}}^{\prime }=\frac{\mathit{\alpha }{Z}_{\mathrm{3}}}{V\mathrm{\Theta }}{e}^{-\frac{\mathrm{1}-X_{\mathrm{3}}}{V\mathrm{\Theta }}},{\mathit{\gamma }}_{\mathrm{2}}^{\prime }=-\frac{\mathit{\alpha }{Z}_{\mathrm{3}}}{V\mathrm{\Theta }}{e}^{-\frac{\mathit{\beta }-Y_{\mathrm{3}}}{V\mathrm{\Theta }}}.$$ Thanks to Eq. (), the coordinates of $P_{\mathrm{3}}$ can be expressed as functions of $Z_{\mathrm{1}}$. The function $T^{\prime }/\mathrm{\Theta }\left(Z_{\mathrm{1}}\right)$ is shown in Fig. 5a. The duration of the interseismic interval ranges from about 2 to 13$\mathrm{\Theta }$, which is from about 60 to 390 a. The maximum value is obtained for $Z_{\mathrm{1}}$ $=$ $Z_{\text{c}}$. We conclude that the evolution of the system after the 1964 event depends on the particular state $P_{\mathrm{1}}$ from which the 1964 event was originated. Since we have expressed $X_{\mathrm{1}}$ and $Y_{\mathrm{1}}$ as functions of $Z_{\mathrm{1}}$, we may characterize the evolution by the value of $Z_{\mathrm{1}}$ as well.

In general, the next event will be an $n$-mode event beginning at a point $P_{\mathrm{4}}$ with coordinates $$X_{\mathrm{4}}=X_{\mathrm{3}}+V{T}^{\prime },Y_{\mathrm{4}}=Y_{\mathrm{3}}+V{T}^{\prime },Z_{\mathrm{4}}=Z_{\mathrm{3}}{e}^{-T^{\prime }/\mathrm{\Theta }},$$ where $T^{\prime }$ is given by Eq. (). There is an infinite number of possible points $P_{\mathrm{4}}$ belonging in part to face $ACD$, in part to $BCD$. Thanks to Eqs. (), () and (), the forces at $P_{\mathrm{4}}$ are $$\begin{array}{cc}{F}_{\mathrm{1}}\left(T_{\mathrm{4}}\right)& =-\mathit{\alpha }{Z}_{\mathrm{1}}-\mathrm{1}+U-V{T}^{\prime }\\ & & +\mathit{\alpha }\left[Z_{\mathrm{1}},+,(,\mathrm{1},-,\mathit{\beta },),U\right]e^{-T^{\prime }/\mathrm{\Theta }},\end{array}$$ $$\begin{array}{cc}{F}_{\mathrm{2}}\left(T_{\mathrm{4}}\right)& =\mathit{\alpha }{Z}_{\mathrm{1}}-\mathit{\beta }+(\mathit{\alpha }+\mathit{\beta })U-V{T}^{\prime }\\ & & -\mathit{\alpha }\left[Z_{\mathrm{1}},+,(,\mathrm{1},-,\mathit{\beta },),U\right]e^{-T^{\prime }/\mathrm{\Theta }}.\end{array}$$ In contrast with the forces in Eq. () at $P_{\mathrm{3}}$, they depend on the particular point $P_{\mathrm{4}}$, hence on $Z_{\mathrm{1}}$ (Fig. 5b), a consequence of viscoelastic relaxation during the interseismic interval.

Hence, the interval [$Z_{\text{a}}$, $Z_{\text{b}}$] can be divided into subintervals leading to different evolutions. If $-$4.55 $\le$ $Z_{\mathrm{1}}$ $<$ 0.20, the next earthquake will be a one-mode event 01. If 0.20 $\le$ $Z_{\mathrm{1}}$ $<$ 0.41, it will be a three-mode event 01–11–10. If $Z_{\mathrm{1}}$ $=$ 0.41, it will be a two-mode event 11–10. If 0.41 $<$ $Z_{\mathrm{1}}$ $<$ 0.70, it will be a three-mode event 10–11–10. Finally, if 0.70 $\le$ $Z_{\mathrm{1}}$ $\le$ 2.86, it will be a one-mode event 10.

The corresponding values of the seismic moment $M_{\mathrm{0}}$ calculated from Eq. () are shown in Fig. 5c and compared with the moment of the 1964 earthquake. It can be seen that the occurrence of an event with a moment greater than the 1964 one is possible only if the value of $Z_{\mathrm{1}}$ is in a narrow range, entailing a narrow range of possible stress distributions on the fault.

Examples of moment rates $\dot{M}$ for the next great Alaska earthquake are shown in Fig. 6 for different values of $Z_{\mathrm{1}}$. The graphs show moment rates for one-mode events 01 (Fig. 6a) and 10 (Fig. 6e), for a two-mode event 11–10 (Fig. 6c), and for three-mode events 01–11–10 (Fig. 6b) and 10–11–10 (Fig. 6d).

Conclusions

We considered a fault with two asperities of different strengths placed in a shear zone subject to a constant strain rate by the motion of adjacent tectonic plates. The equations of motion were written under the hypothesis that the asperities have the same area: this is a reasonable approximation for many earthquakes. The system has been represented by a discrete model described by three variables: the slip deficits of the asperities and the viscoelastic deformation. The system dynamics has one sticking mode and three slipping modes, for which we solved analytically the equations of motion.

If the state of the fault at a given instant in time is known in terms of the system variables, we can calculate the orbit of the system in the phase space and hence predict its evolution. The state of a fault is not directly measurable, but the model shows that the knowledge of the earthquake source functions allows us to constrain the orbit of the system.

The study of the sticking region of the phase space shows how the state of the system before a seismic event controls the sequence of slipping modes in the event. Since the moment rate depends on the number and the sequence of slipping modes, the knowledge of the source function of an earthquake constrains the possible states of the system, hence its orbit in the phase space. Then, if we knew the source functions of a sufficiently large number of consecutive earthquakes, we could constrain the orbit more and more and predict its evolution with a smaller uncertainty.

As an example, we considered the fault that originated the 1964 Alaska earthquake. This earthquake was due to the failure of two distinct asperities; being a large-size event, it was followed by remarkable post-seismic deformation; in addition, more than 50 years have elapsed since the earthquake, allowing such a deformation to be observed over a sufficiently long period of time. The knowledge of the source function of this earthquake allows us to determine the subset of phase space in which the system was before 1964 and the subset to which it came afterwards. This constrains the evolution of the system to a certain bundle of orbits in the phase space but still allows for a wide range of possible occurrence times and source functions for the next earthquake. However, when the next earthquake occurs, the knowledge of its moment rate will allow us to further constrain the orbit, and so on.

The present model is of course a simplification of a real fault, but it suggests how the accumulation of knowledge on the seismic activity of a fault may allow us to constrain the state of the fault and to predict its future activity.

Acknowledgements

The authors are grateful to the editor Ilya Zaliapin, to J. Freymueller and to an anonymous referee for useful comments on the first version of the paper. Edited by: I. Zaliapin Reviewed by: J. Freymueller and one anonymous referee