Methodology

This chapter presents the theory behind PyNetDesign. It describes how the minimum detectable moment magnitude of a microseismic event is evaluated for a network of receivers in a homogeneous, isotropic and attenuating medium, separately for each seismic phase and for each receiver type.

Note

PyNetDesign models a homogeneous medium. The ray path between a source and a receiver is therefore a straight line, there are no interfaces and no transmission losses, and geometrical spreading is simply proportional to distance.

---

Spectral conventions and source-model assumption

We work with frequency in Hz, \(f\), and angular frequency \(\omega = 2\pi f\).

Let \(|U(f)|\) denote the amplitude spectrum magnitude of particle displacement (projected onto the considered phase polarization) recorded at the receiver. If the instrument measures particle velocity, the corresponding spectrum is

(1)\[|V(f)| = 2\pi f\,|U(f)|,\]

which follows from \(v(t)=\dot{u}(t)\).

For shear-dominated microseismic events we adopt the standard \(\omega^{-2}\) source model [Aki and Richards, 2002]. The source displacement spectrum is

\[\begin{split}\Omega_D(f)= \begin{cases} C_0 M_0, & f < f_c, \\ C_0 M_0 \left(\dfrac{f_c}{f}\right)^2, & f \ge f_c, \end{cases}\end{split}\]

where \(M_0\) is seismic moment, \(f_c\) is corner frequency, and \(C_0\) contains radiation-pattern and medium constants.

We define

\[\Omega_0 \equiv C_0 M_0,\]

which is the low-frequency (plateau) level of the far-field source displacement spectrum.

We assume a representative frequency \(f\) below the corner frequency (\(f < f_c\), often \(f=f_{\mathrm{peak}}\), see Section Peak frequency), so the source spectrum is approximately flat, \(\Omega_D(f) \approx \Omega_0\). In this regime, the source spectral function reduces to its low-frequency limit and can be approximated as unity.

---

Intrinsic attenuation

The attenuation operator \(t^*\) (known as global absorption factor [Červený, 2001]) is the path integral along the ray:

(2)\[t^* = \int_R \frac{dr}{v(r)\,Q(r)},\]

For a homogeneous medium, where \(v\) and \(Q\) are constant along the straight ray of length \(r\), Equation (2) reduces to

(3)\[t^* = \frac{r}{vQ}.\]

---

Peak frequency

The peak frequency of an observed signal is the frequency corresponding to the maximum amplitude in Fourier spectra of the particle velocity. As shown by Eisner et al. [2013], the peak frequency can be determined for geophone recordings (i.e., particle velocity) from the global absorption factor \(t^*\) only. The particle velocity is convenient mainly because most of the observations are carried out on instruments recording particle velocity (e.g., geophones or even broadband seismometers).

In Section Spectral conventions and source-model assumption we have shown that the Fourier spectrum of particle displacement is flat at low frequencies with nonzero permanent displacement at the source [Aki and Richards, 2002]. At higher frequencies, the source spectrum is limited by the rise and rupture times of an earthquake fault rupture. These limitations result in a sharp slope change of spectral power density of the particle displacement spectrum. This abrupt change in slope occurs at a corner frequency \(f_c\). Hence, the peak of this positive function must occur for frequencies at or below the corner frequency [Eisner et al., 2013].

The peak frequency \(f_{\mathrm{peak}}\) is the frequency of the maximum amplitude in the particle velocity spectra. Propagation attenuation modifies the observed spectrum through the global absorption factor \(t^*\) (e.g., Červený [2001], Eisner et al. [2013]):

\[|V(f)| \propto f\, e^{-\pi f t^*}.\]

Maximization of this expression for \(f < f_c\) yields

(4)\[f_{\mathrm{peak}} = \frac{1}{\pi t^*},\]

where \(t^*\) is evaluated according to Equation (2).

Note

Equation (4) is only valid for \(f < f_c\), i.e., when the source spectrum is approximately flat. If not, the maximum occurs at the boundary \(f_{\mathrm{peak}}=f_c\).

Similarly, the peak amplitude for particle acceleration will occur at a frequency twice bigger than those for particle velocity, i.e., at \(f_{\mathrm{peak}}=2/(\pi t^*)\).

Taking into account Equation (3), in homogeneous medium the peak amplitude for particle velocity is

(5)\[f_{\mathrm{peak}} = \frac{vQ}{\pi r},\]

where \(v\) is the body wave velocity (P- or S-wave), \(Q\) is the attenuation factor (P- or S-wave), and \(r\) is the distance from a source to a receiver.

---

Detection threshold

An event is considered detected on a receiver when the amplitude of the arriving phase exceeds the local noise level by a required factor. For receiver \(i\) the minimum detectable displacement amplitude spectrum is

(6)\[|U_{\min,i}(f)| = D_{\mathrm{SNR}}\,A_{N,i}\,S(f),\]

where \(A_{N,i}\) is the measured noise level of that receiver, \(D_{\mathrm{SNR}}\) is the signal-to-noise ratio required for detection, and \(S(f)\) converts the recorded quantity to a displacement amplitude spectrum.

Following Equation (1), the conversion coefficient is

Recorded quantity	\(S(f)\)	Comment
displacement	\(1\)	already a displacement
particle velocity	\((2\pi f)^{-1}\)	from \(|V(f)| = 2\pi f\,|U(f)|\)
particle acceleration	\((2\pi f)^{-2}\)	two integrations
strain	\(L_g\)	strain averaged over the gauge length \(L_g\)
strain rate	\(L_g / (2\pi f)\)	as strain, plus one integration in time

For DAS, the recorded quantity is the strain or strain rate along the cable, averaged over the gauge length \(L_g\). Multiplying by \(L_g\) converts it to an equivalent displacement, which assumes the strain to be constant over the gauge length.

The threshold is then corrected for the receiver response of Section Receiver projection and DAS directionality, giving the polarization-aligned threshold \(\widetilde{U}_{\min,i}^{\phi} = C_\phi\,|U_{\min,i}(f)|\) for phase \(\phi\), and divided by the free-surface factor of Section Free surface boundary condition for receivers at the surface.

Network detectability

Inserting the threshold into Equation (7) and converting to a magnitude yields a per-receiver threshold magnitude \(M_w^{\min}(\mathbf{x}, i)\) for every grid point \(\mathbf{x}\) and every receiver \(i\).

An event is only useful if it is seen on several receivers. Sorting the per-receiver values at a grid point in ascending order, the network threshold for a phase is the \(N\)-th smallest value, where \(N\) is the minimum number of receivers required for detection. For joint P- and S-wave detection, the stricter of the two phase thresholds applies,

\[M_w^{\min} = \max\!\left(M_w^{\min,P},\; M_w^{\min,S}\right).\]

For the composite S mode, the SV and SH branches are evaluated separately for each source-receiver pair and the more detectable branch, that is the smaller magnitude, is retained. A phase whose projection vanishes is excluded rather than treated as infinitely detectable.

---

Free surface boundary condition

At the free surface boundary, incoming body waves are reflected, resulting in a complex interference pattern that modifies the recorded amplitude. For the purposes of detectability, we approximate the recorded displacement amplitude to be twice that of the incident wave. This approximation corresponds exactly to the free-surface amplification factor for normal or near-normal incidence, as derived from the free-surface reflection equations (Aki and Richards [2002], Section 5.2.2).

---

Receiver projection and DAS directionality

PyNetDesign treats ordinary geophones as three-component receivers by default. For these receivers the recorded scalar amplitude is assumed to be aligned with the arriving phase polarization, so no projection loss is applied. If a station geometry row has Components=Z, the arriving displacement is projected onto the vertical axis. For DAS receivers, PyNetDesign applies the same phase-dependent projection logic onto the local fiber tangent \(\boldsymbol\tau\).

In a homogeneous medium the ray is a straight line, so the phase polarization is computed from the straight receiver-to-grid direction. The amplitude correction for phase \(\phi\) and receiver direction \(\mathbf d\) is

\[C_\phi = \frac{1}{|\mathbf d \cdot \mathbf e_\phi|},\]

where \(\mathbf d\) is \(\hat{\mathbf z}\) for vertical-component geophones and \(\boldsymbol\tau\) for DAS. \(\mathbf e_P\) is the longitudinal P polarization, \(\mathbf e_{SV}\) is the transverse polarization in the local vertical ray plane, and \(\mathbf e_{SH}\) is perpendicular to that plane. The corrected minimum recorded amplitude is multiplied by \(C_\phi\). For composite S mode, SV and SH are corrected separately and the more detectable branch is retained.

If the projection is zero, that phase is not observable on the selected receiver component for that source-station pair. For example, SH is horizontal in the current isotropic convention and is therefore not recorded by a vertical-component geophone.

For nearly vertical rays, the explicit SV/SH azimuthal basis is convention-dependent because the horizontal ray azimuth is poorly defined. PyNetDesign uses a stable horizontal fallback basis in this case; the composite S mode is the practical robust choice when only S-wave detectability is required.

---

Radiation pattern

In the far field, body-wave amplitudes from a seismic moment tensor source depend on the take-off direction and on the wave polarization. Let \(\mathbf n\) be the unit vector in the P-wave take-off direction at the source (pointing from source to receiver), and let \(\mathbf e\) be the unit polarization vector of the considered phase at the source.

Following the quantitative formulation of Aki and Richards [2002] (Section 4.5), for a moment tensor \(\mathbf M\), the dimensionless scalar radiation pattern factor for the phase with polarization \(\mathbf e\) is

\[R(\mathbf n,\mathbf e) = \frac{\mathbf e^{T}\,\mathbf M\,\mathbf n}{M_0},\]

where \(M_0\) is the scalar seismic moment (used only for normalization of the dimensionless factor).

For P waves, the polarization is longitudinal (\(\mathbf e=\mathbf n\)), leading to the familiar radiation pattern factor (Aki and Richards [2002], Eq. 4.89):

\[R_P = \frac{\mathbf n^{T}\,\mathbf M\,\mathbf n}{M_0}.\]

For S waves, there are two orthogonal transverse polarizations. Denoting by \(\mathbf s_{SV}\) the unit SV polarization (in the ray plane) and by \(\mathbf s_{SH}\) the unit SH polarization (perpendicular to the ray plane), the corresponding radiation factors are (Aki and Richards [2002], Eqs. 4.90 and 4.91):

\[R_{SV} = \frac{\mathbf s_{SV}^{T}\,\mathbf M\,\mathbf n}{M_0}, \qquad R_{SH} = \frac{\mathbf s_{SH}^{T}\,\mathbf M\,\mathbf n}{M_0}.\]

In this documentation and in the magnitude-sensitivity formulas, \(|R|\) denotes the magnitude of the appropriate factor (\(|R_P|\) for P, \(|R_{SV}|\) for SV, and \(|R_{SH}|\) for SH). For wave_mode="S", PyNetDesign evaluates SV and SH separately and keeps the more detectable branch for each source-station pair.

For microseismic network design, the source mechanism and take-off directions are often unknown. In that case, it is common to use an effective radiation-pattern magnitude (i.e., a somehow averaged value of \(|R|\) over a suitable ensemble of mechanisms and directions). Hallo and Eisner [2013] following Boore and Boatwright [1984] suggest to use representative RMS values for DC seismic sources \(R_P^{\mathrm{eff}} = \sqrt{\langle R_P^2\rangle} \approx 0.52\) and \(R_S^{\mathrm{eff}} = \sqrt{\langle R_S^2\rangle} \approx 0.63\). PyNetDesign uses fixed RMS constants \(R_P^{\mathrm{eff}}=\sqrt{4/15}\approx 0.516\), \(R_{SV}^{\mathrm{eff}}=\sqrt{7/30}\approx 0.483\) for SV and \(R_{SH}^{\mathrm{eff}} = \sqrt{1/6}\approx 0.408\) for SH. See Section Derivation of effective radiation pattern magnitudes for more information.

---

Magnitude sensitivity

The magnitude sensitivity module is designed to estimate the magnitudes of microseismic events that can be theoretically detected by the network of receivers. This sensitivity (or detectability) varies in space and hence it is presented in the form of horizontal or vertical slices in the 3D space. PyNetDesign uses attenuation of the dominant signal (based on the quality factor of the medium \(Q\) ), and signal-to-noise ratio (SNR) necessary for seismic event detection. Results depend on a model of attenuation, noise, wave velocity, and detection method. The resulting sensitivity can show the spatial distribution of the minimal technically detectable magnitude, or, in other words, the spatial distribution of the magnitude of detection catalogue completeness. The moment magnitude of an event, according to Kanamori [1977], is given by logarithm of the scalar value of the seismic moment:

\[M_w = \frac{2}{3}\left(\log_{10}(M_0)-9.1\right),\]

where \(M_0\) is the seismic moment.

PyNetDesign computes minimal moment magnitude of an event, which will theoretically produce amplitudes on receivers necessary for the event detection. Amplitude levels for the event detection are dependent on local noise levels and signal-to-noise ratio necessary for particular detection method.

Homogeneous medium

In a homogeneous isotropic medium, the scalar seismic moment \(M_0\) can be related to the magnitude of the displacement Fourier spectrum \(|U(f)|\) (projected onto the arriving polarization of the considered phase P or S) at frequency \(f\) as

(7)\[M_0 = \frac{4 \pi \rho v^3 r}{|R|} \;|U(f)|\; \exp\!\big(\pi f\,t^*\big),\]

where \(\rho\) is density, \(v\) is the body-wave velocity (\(v=\alpha\) for P waves, \(v=\beta\) for S waves), and \(r\) is the source-receiver distance. The factor \(|R|\) is the dimensionless radiation-pattern magnitude for the considered phase, consistent with the moment-tensor radiation convention of Aki and Richards [2002] (see Section Radiation pattern).

For derivation of formula (7), see Section Derivation of the seismic moment in a homogeneous medium.

The exponential term \(\exp(\pi f t^*)\) corrects for intrinsic attenuation using the attenuation operator \(t^*\) (see Section Intrinsic attenuation).

Because P-wave velocity is typically higher than S-wave velocity, for the same propagation path and \(Q\) the term \(t^*=r/(vQ)\) is smaller for P waves, so the attenuation correction \(\exp(\pi f t^*)\) is closer to unity. In practice, detectability is also controlled by radiation pattern and noise characteristics; network-design sensitivity is therefore often evaluated separately for P and S phases [Hallo and Eisner, 2013].

---

Derivation of the seismic moment in a homogeneous medium

Here we derive Equation (7) from the far-field radiation of a point moment-tensor source.

In an unbounded homogeneous isotropic medium, the far-field displacement produced by a point source with moment time function \(M_0(t)\) is, for the phase with radiation-pattern factor \(R\) and velocity \(v\) (Aki and Richards [2002], Section 4.3),

(8)\[u(t) = \frac{R}{4\pi\rho v^3 r}\,\dot{M}_0\!\left(t - \frac{r}{v}\right),\]

where the \(1/r\) factor is the geometrical spreading of a spherical wave in a homogeneous medium and the dot denotes the time derivative.

Taking the Fourier transform and using the fact that the moment rate spectrum tends to the scalar seismic moment at low frequency, \(|\dot{M}_0(f)| \to M_0\) for \(f \ll f_c\), the low-frequency plateau of the displacement amplitude spectrum is

(9)\[\Omega_0 = \frac{M_0\,|R|}{4\pi\rho v^3 r},\]

consistent with the source-model assumption of Section Spectral conventions and source-model assumption.

Intrinsic attenuation multiplies the spectrum by the dissipation filter of Section Intrinsic attenuation, so the spectrum actually observed at the receiver is

(10)\[|U(f)| = \frac{M_0\,|R|}{4\pi\rho v^3 r}\,\exp\!\big(-\pi f\,t^*\big), \qquad t^* = \frac{r}{vQ}.\]

Solving Equation (10) for the seismic moment gives

\[M_0 = \frac{4\pi\rho v^3 r}{|R|}\; |U(f)|\; \exp\!\left(\pi f\,\frac{r}{vQ}\right),\]

which is Equation (7).

Note

The exponential appears with a positive exponent in the moment formula: the amplitude observed at the receiver has already been reduced by attenuation, so recovering the source moment means undoing that reduction.

Note also that \(|U(f)|\) is an amplitude spectrum. The conversion of a recorded noise level to \(|U(f)|\) is carried out once, by the coefficient \(S(f)\) of Section Detection threshold, and must not be applied a second time here.

---

Derivation of effective radiation pattern magnitudes

For microseismic network design with unknown focal mechanisms, we utilize effective averaged radiation factors. While literature often name these as average magnitudes, they are mathematically evaluated as the Root Mean Square (RMS) averages over the entire focal sphere to preserve wave energy [Boore and Boatwright, 1984].

We assume a canonical double-couple (DC) source. In a coordinate system where the fault plane normal is aligned with the x-axis and the slip direction is aligned with the y-axis, the seismic moment tensor \(\mathbf{M}\) is:

\[\begin{split}\mathbf{M} = M_0 \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}.\end{split}\]

Let \(\mathbf{\gamma}\) be the unit take-off vector parameterized by polar angle \(\theta\) and azimuthal angle \(\phi\):

\[\begin{split}\mathbf{\gamma} = \begin{pmatrix} \sin\theta\cos\phi \\ \sin\theta\sin\phi \\ \cos\theta \end{pmatrix}.\end{split}\]

P-wave average

The P-wave radiation pattern factor is the projection of the equivalent force \(\mathbf{M}\mathbf{\gamma}\) onto the longitudinal polarization \(\mathbf{\gamma}\):

\[R_P = \frac{\mathbf{\gamma}^T \mathbf{M} \mathbf{\gamma}}{M_0} = 2 \sin^2\theta \sin\phi \cos\phi = \sin^2\theta \sin 2\phi.\]

To find the RMS average, we integrate the square of \(R_P\) over the unit sphere (where the averaging operator is \(\frac{1}{4\pi} \int_0^{2\pi} \int_0^\pi \dots \sin\theta \, d\theta \, d\phi\)):

\[\langle R_P^2 \rangle = \frac{1}{4\pi} \int_0^{2\pi} \int_0^\pi \left(\sin^2\theta \sin 2\phi\right)^2 \sin\theta \, d\theta \, d\phi.\]

Separating the integrals:

\[\langle R_P^2 \rangle = \frac{1}{4\pi} \left( \int_0^\pi \sin^5\theta \, d\theta \right) \left( \int_0^{2\pi} \sin^2 2\phi \, d\phi \right).\]

Evaluating the integrals yields \(16/15\) for the polar part and \(\pi\) for the azimuthal part:

\[\langle R_P^2 \rangle = \frac{1}{4\pi} \left(\frac{16}{15}\right) (\pi) = \frac{4}{15}.\]

Taking the square root provides the effective P-wave magnitude:

\[R_P^{\mathrm{rms}} = \sqrt{\frac{4}{15}} \approx 0.516 \approx 0.52.\]

S-wave average

For S-waves, the total radiation pattern squared, \(R_S^2\), is the squared magnitude of the transverse component of the equivalent force. This is the total equivalent force squared minus the longitudinal (P-wave) component squared:

\[R_S^2 = \frac{|\mathbf{M}\mathbf{\gamma}|^2}{M_0^2} - R_P^2.\]

First, we compute the squared magnitude of the normalized equivalent force:

\[\frac{|\mathbf{M}\mathbf{\gamma}|^2}{M_0^2} = (\sin\theta\sin\phi)^2 + (\sin\theta\cos\phi)^2 + 0^2 = \sin^2\theta.\]

Averaging this over the focal sphere:

\[\langle \frac{|\mathbf{M}\mathbf{\gamma}|^2}{M_0^2} \rangle = \frac{1}{4\pi} \int_0^{2\pi} \int_0^\pi \sin^2\theta \sin\theta \, d\theta \, d\phi = \frac{1}{4\pi} \left( \frac{4}{3} \right) (2\pi) = \frac{2}{3}.\]

Subtracting the P-wave averaged from the total force average yields the total S-wave average:

\[\langle R_S^2 \rangle = \frac{2}{3} - \langle R_P^2 \rangle = \frac{10}{15} - \frac{4}{15} = \frac{6}{15} = \frac{2}{5}.\]

Taking the square root provides the effective total S-wave magnitude:

\[R_S^{\mathrm{rms}} = \sqrt{\frac{2}{5}} \approx 0.632 \approx 0.63.\]

SV- and SH-wave average

The total S-wave coefficient derived above describes the magnitude of the full transverse radiation vector. It does not imply equal partition between the SV and SH components.

To obtain component-wise averages, we now switch to the standard fault-based angular convention of Aki and Richards [2002]. To avoid confusion with the spherical angles used above in the tensor-based derivation, we denote these fault-based angles by \(\Theta\) and \(\Phi\).

In this convention, the standard double-couple far-field radiation-pattern expressions for P, SV, and SH waves are:

\[R_P(\Theta,\Phi) = \sin 2\Theta \cos\Phi,\]

\[R_{SV}(\Theta,\Phi) = \cos 2\Theta \cos\Phi,\]

\[R_{SH}(\Theta,\Phi) = -\cos\Theta \sin\Phi,\]

where \(\Theta\) is the take-off angle and \(\Phi\) is the azimuth in the fault-based coordinate system.

The effective component-wise coefficients are defined by RMS averaging over the focal sphere:

\[\langle R_{SV}^2 \rangle = \frac{1}{4\pi} \int_0^{2\pi}\int_0^\pi \cos^2 2\Theta \cos^2\Phi \sin\Theta\,d\Theta\,d\Phi,\]

\[\langle R_{SH}^2 \rangle = \frac{1}{4\pi} \int_0^{2\pi}\int_0^\pi \cos^2\Theta \sin^2\Phi \sin\Theta\,d\Theta\,d\Phi.\]

Using

\[\int_0^{2\pi}\cos^2\Phi\,d\Phi = \pi, \qquad \int_0^{2\pi}\sin^2\Phi\,d\Phi = \pi,\]

\[\int_0^\pi \cos^2 2\Theta \sin\Theta\,d\Theta = \frac{14}{15}, \qquad \int_0^\pi \cos^2\Theta \sin\Theta\,d\Theta = \frac{2}{3},\]

we obtain

\[\langle R_{SV}^2 \rangle = \frac{1}{4\pi}\,\pi\,\frac{14}{15} = \frac{7}{30},\]

\[\langle R_{SH}^2 \rangle = \frac{1}{4\pi}\,\pi\,\frac{2}{3} = \frac{1}{6}.\]

Therefore, the RMS effective coefficients are

\[R_{SV}^{\mathrm{rms}} = \sqrt{\langle R_{SV}^2\rangle} = \sqrt{\frac{7}{30}} \approx 0.483,\]

\[R_{SH}^{\mathrm{rms}} = \sqrt{\langle R_{SH}^2\rangle} = \sqrt{\frac{1}{6}} \approx 0.408.\]

As expected, the total S-wave mean-square coefficient is the sum of the two orthogonal transverse components:

\[\langle R_S^2\rangle = \langle R_{SV}^2\rangle + \langle R_{SH}^2\rangle = \frac{7}{30} + \frac{1}{6} = \frac{2}{5}.\]

Hence

\[R_{S}^{\mathrm{rms}} = \sqrt{\langle R_S^2\rangle} = \sqrt{\frac{2}{5}} \approx 0.63.\]

Note

Therefore, the commonly used value \(0.63\) corresponds to the combined S-wave amplitude, not to SV alone and not to SH alone. PyNetDesign uses fixed RMS constants \(R_P^{\mathrm{rms}}=\sqrt{4/15}\), \(R_{SV}^{\mathrm{rms}}=\sqrt{7/30}\) for explicit SV, \(R_{SH}^{\mathrm{rms}}=\sqrt{1/6}\) for explicit SH, and evaluates both S branches for composite S detectability.

---

References

[1] (1,2,3,4,5,6,7,8,9)

K. Aki and P. G. Richards. Quantitative Seismology. University Science Books, 2002. ISBN 9780935702965.

[2] (1,2)

V. Červený. Seismic Ray Theory. Cambridge University Press, 2001. doi:10.1017/CBO9780511529399.

[3] (1,2,3)

L. Eisner, D. Gei, M. Hallo, I. Opršal, and M. Y. Ali. The peak frequency of direct waves for microseismic events. GEOPHYSICS, 78(6):A45–A49, November 2013. doi:10.1190/geo2013-0197.1.

[4] (1,2)

M. Hallo and L. Eisner. Influence of the double-couple source model on the sensitivity of microseismic monitoring networks. In London 2013, 75th eage conference en exhibition incorporating SPE Europec, London2013. EAGE Publications BV, 2013. doi:10.3997/2214-4609.20131023.

[5] (1,2)

David M. Boore and John Boatwright. Average body-wave radiation coefficients. Bulletin of the Seismological Society of America, 74(5):1615–1621, October 1984. doi:10.1785/bssa0740051615.

[6]

H. Kanamori. The energy release in great earthquakes. Journal of Geophysical Research, 82(20):2981–2987, 1977. doi:10.1029/jb082i020p02981.