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 possibly reject low SNR spectra) and, optionally, to weight the spectral inversion.
The component spectra are combined through the root sum of squares:
where \(f\) is the frequency and \(S_x(f)\) is the P- or S-wave spectrum for component \(x\).
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):
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.
Following 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:
with
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:
Then the spectrum is converted in units of magnitude (the \(Y_x (f)\) vector used in the inversion):
The data vector is compared to the theoretical model:
Finally coming to the following model used for the inversion:
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:
TNC: truncated Newton algorithm (with bounds)
LM: Levenberg-Marquardt algorithm (warning: Trust Region Reflective algorithm will be used instead if bounds are provided)
GS: grid search
IS: importance sampling of misfit grid, using k-d tree
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 measured 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 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.
Attenuation
The retrieved attenuation parameter \(t^*\) is converted to the P- or S-wave quality factor \(Q_0^{[P|S]}\) using the following expression:
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).