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 (e.g., Z, N, E), within a predefined time window.

The same thing is done for a noise time window: noise spectrum is used to compute spectral signal-to-noise ratio (and possibily reject low SNR spectra) and, optionnaly, 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:

\[S(f) = \sqrt{S^2_z(f) + S^2_n(f) + S^2_e(f)}\]

where \(f\) is the frequency and \(S_x(f)\) is the P- or S-wave spectrum for component \(x\).

Example spectrum plot

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

Spectral model

The Fourier amplitude spectrum of the 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(f) = \frac{1}{\mathcal{G}(r)} \times \frac{2 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_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;

  • \(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_c\) is the corner frequency;

  • \(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 in one of the following ways:

  • \(\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.

  • Follwing Boatwright et al. [2002], eq. 8:

    • 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 (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.

Building spectra

In source_spec, the observed spectrum of component \(x\), \(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_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_x(f) = M_O \times \frac{1}{1+\left(\frac{f}{f_c}\right)^2} \times e^{- \pi f t^*}\]

Then the spectrum is converted in units of magnitude (the \(Y_x (f)\) vector used in the inversion):

\[Y_x(f) \equiv \frac{2}{3} \times \left( \log_{10} M_x(f) - 9.1 \right)\]

The data vector is compared to the teoretical model:

\[ \begin{align}\begin{aligned}Y_x(f) = \frac{2}{3} \left[ \log_{10} \left( M_O \times \frac{1}{1+\left(\frac{f}{f_c}\right)^2} \times e^{- \pi f t^*} \right) - 9.1 \right] =\\ = \frac{2}{3} (\log_{10} M_0 - 9.1) + \frac{2}{3} \left[ \log_{10} \left( \frac{1}{1+\left(\frac{f}{f_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_x(f) = M_w + \frac{2}{3} \left[ - \log_{10} \left( 1+\left(\frac{f}{f_c}\right)^2 \right) - \pi \, f t^* \log_{10} e \right]\]

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

Inversion procedure

The parameters to determine are \(M_w\), \(f_c\) and \(t^*\).

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

Starting from the inverted parameters \(M_0\) ( \(M_w\) ), \(fc\), \(t^*\) and following the equations in Madariaga [2011], other quantities are computed for each station:

  • the Brune stress drop

  • the source radius

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

Finally, the radiated energy \(E_r\) can be mesured from the displacement spectra, following the approach described in Lancieri et al. [2012].

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

Event averages are computed from single station estimates. Outliers are rejected based on the interquartile range rule.

Attenuation

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-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).