Theoretical Background

Overview

source_spec inverts the P- or S-wave displacement amplitude spectra from station recordings of a single event.

For each station, the code computes P- or S-wave displacement amplitude spectra for each component \(x\) (e.g., Z, N, E), within a predefined time window.

\[S_x^{p|s}(f) = \left| \int_{t^{p|s}_0}^{t^{p|s}_1} s_x(t) e^{-i 2 \pi f t} dt \right|\]

where the exponent \(p|s\) means that we are considering either P- or S-waves, \(t^{p|s}_0\) and \(t^{p|s}_1\) are the start and end times of the P- or S-wave time window, \(s_x(t)\) is the displacement time series for component \(x\), and \(f\) is the frequency.

Note that the Fourier amplitude spectrum of ground displacement has the dimensions of displacement (\(s_x(t)\)) multiplied by time (\(dt\)).

The same thing is done for a noise time window: noise spectrum is used to compute spectral signal-to-noise ratio (and possibly reject low SNR spectra) and, optionally, to weight the spectral inversion.

Example trace plot

Example three-component trace plot (in velocity), showing noise and S-wave windows.

The component spectra are combined through the root sum of squares (e.g., Z, N, E):

\[S^{p|s}(f) = \sqrt{ \left( S^{p|s}_z(f) \right)^2 + \left( S^{p|s}_n(f) \right)^2 + \left( S^{p|s}_e(f) \right)^2 }\]

(This is actually done later in the code, after converting the spectra to magnitude units, see below.)

Example spectrum plot

Example displacement spectrum for noise and S-wave, including inversion results.

Spectral model

The Fourier amplitude spectrum of the P- or S-wave displacement in far field can be modelled as the product of a source term [Brune, 1970] and a propagation term (geometric and anelastic attenuation of body waves):

\[S^{p|s}(f) = \frac{1}{\mathcal{G}(r)} \times \frac{F R_{\Theta\Phi}} {4 \pi \rho_h^{1/2} \rho_r^{1/2} c_h^{5/2} c_r^{1/2}} \times M_O \times \frac{1}{1+\left(\frac{f}{f^{p|s}_c}\right)^2} \times e^{- \pi f t^*}\]

where:

  • \(\mathcal{G}(r)\) is the geometrical spreading coefficient (see below) and \(r\) is the hypocentral distance;

  • \(F\) is the free surface amplification factor (generally assumed to be \(2\));

  • \(R_{\Theta\Phi}\) is the radiation pattern coefficient for P- or S-waves (average or computed from focal mechanism, if available);

  • \(\rho_h\) and \(\rho_r\) are the medium densities at the hypocenter and at the receiver, respectively;

  • \(c_h\) and \(c_r\) are the P- or S-wave velocities at the hypocenter and at the receiver, respectively;

  • \(M_O\) is the seismic moment;

  • \(f\) is the frequency;

  • \(f^{p|s}_c\) is the corner frequency for P- or S-waves;

  • \(t^*\) is an attenuation parameter which includes anelastic path attenuation (quality factor) and station-specific effects.

Geometrical spreading

The geometrical spreading coefficient \(\mathcal{G}(r)\) can be defined, for local and regional distances, in one of the following ways (see the geom_spred_model option in Configuration File):

  • \(\mathcal{G}(r) = r^n\): \(n\) can be any positive number. \(n=1\) (default value) is the theoretical value for a body wave in a homogeneous full-space; \(n=0.5\) is the theoretical value for a surface wave in a homogeneous half-space.

  • Following Boatwright et al. [2002] (eq. 8), to account for the mixture of body waves, Lg waves and surface waves at regional distances (\(r < 200 km\)), a two-part geometrical spreading coefficient:

    • body wave spreading (\(\mathcal{G}(r) = r\)) for hypocentral distances below a cutoff distance \(r_0\);

    • frequency dependent spreading for hypocentral distances above the cutoff distance \(r_0\).

More precisely, the expression derived from Boatwright et al. [2002] is:

\[\begin{split}\mathcal{G}(r) = \begin{cases} r & r \le r_0\\ r_0 (r/r_0)^{\gamma (f)} & r > r_0 \end{cases}\end{split}\]

with

\[\begin{split}\gamma (f) = \begin{cases} 0.5 & f \le 0.20 Hz\\ 0.5 + 2 \log_{10} (5f) & 0.20 < f < 0.25 Hz\\ 0.7 & f \ge 0.25 Hz\\ \end{cases}\end{split}\]

Note that here we use the square root of eq. 8 in Boatwright et al. [2002], since we correct the spectral amplitude and not the energy.

For teleseismic distances (see the option geom_spread_min_teleseismic_distance in Configuration File), the geometrical spreading coefficient is defined as in Okal [1992] (eq. 4):

\[\mathcal{G}(\Delta) = \frac{a}{g(\Delta)}\]

where \(\Delta\) is the great circle distance between the source and the receiver, \(a\) is the Earth radius and \(g(\Delta)\) is defined as:

\[g(\Delta) = \left( \frac{\rho_h c_h}{\rho_r c_r} \frac{\sin i_h}{\sin \Delta} \frac{1}{\cos i_r} \left| \frac{d i_h}{d \Delta} \right| \right)^{1/2}\]

where \(\rho_h\) and \(\rho_r\) are the medium densities at the hypocenter and at the receiver, respectively, \(c_h\) and \(c_r\) are the P- or S-wave velocities at the hypocenter and at the receiver, respectively, \(i_h\) and \(i_r\) are the takeoff angle (hypocenter) and the incidence angle (receiver), respectively, and \(\frac{d i_h}{d \Delta}\) is the variation of the takeoff angle within a ray tube of width \(\Delta\) (see Okal [1992] for details).

Building spectra

In source_spec, the observed spectrum of component \(x\) (vertical or horizontal), \(S^{p|s}_x(f)\) is converted into moment magnitude units \(M_w\).

The first step is to multiply the spectrum for the geometrical spreading coefficient and convert it to seismic moment units:

\[M^{p|s}_x(f) \equiv \mathcal{G}(r) \times \frac{4 \pi \rho_h^{1/2} \rho_r^{1/2} c_h^{5/2} c_r^{1/2}} {2 R_{\Theta\Phi}} \times S^{p|s}_x(f) = M_O \times \frac{1}{1+\left(\frac{f}{f^{p|s}_c}\right)^2} \times e^{- \pi f t^*}\]

Then the spectrum is converted in units of magnitude:

\[Y^{p|s}_x(f) \equiv \frac{2}{3} \times \left( \log_{10} M^{p|s}_x(f) - 9.1 \right)\]

And the final data vector \(Y^{p|s}(f)\) is obtained by combining the three components (e.g., Z, N, E) through the root sum of squares:

\[Y^{p|s}(f) = \sqrt{ \left( Y^{p|s}_z(f) \right)^2 + \left( Y^{p|s}_n(f) \right)^2 + \left( Y^{p|s}_e(f) \right)^2 }\]

The data vector is compared to the theoretical model:

\[ \begin{align}\begin{aligned}Y^{p|s}(f) = \frac{2}{3} \left[ \log_{10} \left( M_O \times \frac{1}{1+\left(\frac{f}{f^{p|s}_c}\right)^2} \times e^{- \pi f t^*} \right) - 9.1 \right] =\\ = \frac{2}{3} \left( \log_{10} M_0 - 9.1 \right) + \frac{2}{3} \left[ \log_{10} \left( \frac{1}{1+\left(\frac{f}{f^{p|s}_c}\right)^2} \right) + \log_{10} \left( e^{- \pi f t^*} \right) \right]\end{aligned}\end{align} \]

Finally coming to the following model used for the inversion:

\[Y^{p|s}(f) = M_w + \frac{2}{3} \left[ - \log_{10} \left( 1+\left(\frac{f}{f^{p|s}_c}\right)^2 \right) - \pi \, f t^* \log_{10} e \right]\]

where \(M_w \equiv \frac{2}{3} (\log_{10} M_0 - 9.1)\).

Inverted parameters

The parameters determined from the spectral inversion are \(M_w\), \(f^{p|s}_c\) and \(t^*\).

The inversion is performed in moment magnitude \(M_w\) units (logarithmic amplitude). Different inversion algorithms can be used:

Other computed parameters

Starting from the inverted parameters \(M_0\) ( \(M_w\) ), \(f^{p|s}_c\), \(t^*\) and following the equations in Madariaga [2011] and Lancieri et al. [2012], other quantities are computed for each station:

  • the static stress drop \(\Delta \sigma\)

  • the source radius \(a\)

  • the radiated energy \(E_r\)

  • the apparent stress \(\sigma_a\)

  • the quality factor \(Q_0\) of P- or S-waves

As a bonus, local magnitude \(M_l\) can be computed as well.

Event summaries (mean, weighted mean, percentiles) are computed from single station estimates. For mean and weighted mean estimation, outliers are rejected based on the interquartile range rule.

Source radius

The source radius is computed, assuming a circular rupture model, from the corner frequency \(f^{p|s}_c\) (Kaneko and Shearer [2014], equation 2):

\[a = k^{p|s} \frac{\beta_h}{f^{p|s}_c}\]

where \(\beta_h\) is the shear wave speed at the hypocenter (in \(m/s\)), \(f^{p|s}_c\) is the corner frequency (in \(Hz\)) estimated from the spectral inversion of P or S waves and \(k^{p|s}\) is a constant which depends on the source model.

Brune [1970] provides an expression for \(k^s\) in the case of a static circular crack (equation 31 in Madariaga [2011]):

\[k^s_{Brune} = 0.3724\]

Kaneko and Shearer [2014] compiled a table including their own values for \(k^p\) and \(k^s\) as well as values obtained from other authors. The values are given as a function of the rupture velocity \(V_r\) of a dynamic circular crack (or a static crack, when \(V_r\) is infinite):

Table 1 of Kaneko and Shearer [2014]

\(V_r/\beta_h\)

\(k^p_{K\&S}\)

\(k^s_{K\&S}\)

\(k^p_{Mada}\)

\(k^s_{Mada}\)

\(k^s_{Brune}\)

\(k^p_{S\&H}\)

\(k^s_{S\&H}\)

Infinite

0.3724

0.9

0.38

0.26

0.32

0.21

0.42

0.29

0.8

0.35

0.26

0.39

0.28

0.7

0.32

0.26

0.36

0.27

0.6

0.30

0.25

0.34

0.27

0.5

0.28

0.22

0.31

0.24

Where “K&S” stands for Kaneko and Shearer [2014], “Mada” for Madariaga [1976], and “S&H” for Sato and Hirasawa [1973].

Static stress drop

The static stress drop \(\Delta \sigma\) is computed under the assumption of a circular rupture of radius \(a\), as discussed in Madariaga [2011] (equation 27):

\[\Delta \sigma = \frac{7}{16} \frac{M_0}{a^3}\]

where \(M_0\) is the seismic moment (in \(N \cdot m\)) and \(a\) is the source radius (in \(m\)).

Radiated energy

The computation of the radiated energy \(E_r\) starts with the integral of the squared velocity amplitude spectrum: \([\dot{S}^{p|s}(f)]^2 = [ 2 \pi f S^{p|s}(f) ]^2\).

Following Boatwright et al. [2002] (equation 1) and Lancieri et al. [2012] (equation 3), the P- or S-wave radiated energy is computed as:

\[\tilde{\tilde{E}}_r^{p|s} = 8 \pi \mathcal{G}^2(r) C^2 \rho_r c_r \int_{f_{min}}^{f_{max}} e^{2 \pi f t^*} [\dot{S}^{p|s}(f)]^2 df\]

where \(\mathcal{G}^2(r)\) is the squared geometrical spreading coefficient (see above), \(C\) is a constant discussed below, \(\rho_r\) and \(c_r\) are, respectively, the density and P- or S-wave velocity at the receiver (their product is the seismic impedance), \(f_{min}\) and \(f_{max}\) are the minimum and maximum frequency used to compute the energy (see Configuration File for details on the Er_freq_range parameter), and the exponential term in the integrand is the squared correction for anelastic attenuation. The double tilde on top of \(\tilde{\tilde{E}}_r^{p|s}\) means that the radiated energy needs to be further corrected for noise and finite bandwidth (see below).

The constant \(C\) is defined in Boatwright et al. [2002] (equation 2) as:

\[C = \frac{\left<R_{\Theta\Phi}\right>}{R_{\Theta\Phi} F}\]

where \(\left<R_{\Theta\Phi}\right>\) is the root mean square P- or S-wave radiation pattern computed on the focal sphere, \(R_{\Theta\Phi}\) is the radiation pattern coefficient for the given station, and \(F\) is the free surface amplification factor. If a focal mechanism is not available, then it is assumed \(R_{\Theta\Phi} = \left<R_{\Theta\Phi}\right>\) and, hence, \(C = 1/F\). This assumption means that we rely on the averaging between measurements of radiated energy at different stations, instead of precise measurements at a single station.

Noise correction

To account for low frequency noise, below the corner frequency, under the hypothesis that energy is additive and that noise is stationary, we compute a noise-corrected energy as:

\[\tilde{E}^{p|s}_r = \tilde{\tilde{E}}^{p|s}_r - \tilde{\tilde{E}}^{noise}_r\]

where the first term is the radiated energy computed from the P- or S-wave spectrum and the second term is the radiated energy computed from the noise spectrum. If the above difference is negative, then the measure is rejected, since the noise is too large compared to the signal.

Finite bandwidth correction

Next step is to correct the radiated energy for the missing energy above \(f_{max}\), not taken into account in the integral of the squared velocity amplitude spectrum (finite bandwidth correction). Following Lancieri et al. [2012] (equation 4), and Di Bona and Rovelli [1988], the noise-corrected radiated energy is divided by the theoretical ratio \(R\) between the estimated radiated energy and the true radiated energy, defined as:

\[R = \frac{2}{\pi} \left[ \arctan(f_{max}/f^{p|s}_c) - \frac{f_{max}/f^{p|s}_c}{1+(f_{max}/f^{p|s}_c)^2} \right]\]

where \(f_{max}\) is the maximum frequency used to compute the energy integral and \(f^{p|s}_c\) is the P- or S-wave corner frequency.

The values of R range between 0 (for \(f_{max}/f^{p|s}_c \to 0\)) and 1 (for \(f_{max}/f^{p|s}_c \to \infty\)).

The corrected radiated energy for P- or S-waves is then:

\[E^{p|s}_r = \frac{\tilde{E}^{p|s}_r}{R}\]

Energy partition

The final step is to account for the partition of energy between P and S waves. Following Boatwright and Choy [1986] (equations 8 and 15) the ratio between the radiated energy measured from S-waves and the radiated energy measured from P-waves is:

\[\frac{E_r^s}{E_r^p} = 15.6\]

The final estimate of the radiated energy is then:

\[E_r = \left( 1 + 15.6 \right) E_r^p\]

or

\[E_r = \left( 1 + \frac{1}{15.6} \right) E_r^s\]

depending on whether the radiated energy is computed from P or S waves.

Apparent stress

The apparent stress \(\sigma_a\) is computed as (Madariaga [2011], eq. 18):

\[\sigma_a = \mu_h \frac{E_r}{M_0}\]

where \(\mu_h\) is the shear modulus (or rigidity, in \(Pa\)) near the hypocenter, \(E_r\) is the radiated energy (in \(N \cdot m\)), and \(M_0\) is the seismic moment (in \(N \cdot m\)).

The value of \(\mu_h\) is computed from the shear wave velocity (\(\beta_h\)) and the density (\(\rho_h\)) at the hypocenter, using the following expression:

\[\mu_h = \rho_h \beta_h^2\]

Quality factor

The retrieved attenuation parameter \(t^*\) is converted to the P- or S-wave quality factor \(Q_0^{p|s}\) using the following expression:

\[Q_0^{p|s} = \frac{tt^{p|s}(r)}{t^*}\]

where \(tt^{p|s}(r)\) is the P- or S-wave travel time from source to station and \(r\) is the hypocentral distance.

Station Residuals

Station-specific effects can be determined by running source_spec on several events and computing the average of station residuals between observed and inverted spectra. These averages are obtained through the command source_residuals; the resulting residuals file can be used for a second run of source_spec (see the residuals_filepath option in Configuration File).