Earthquake Magnitudes: Seismic Moment, Moment Magnitude and Local Magnitude

Introduction

There exist several magnitude scales that measure the size of an earthquake. The local magnitude \(M_L\), or Richter’s magnitude, is often used for small events because it is easy to calculate (however, uncertainties are often large!). The moment magnitude \(M_w\) is the only magnitude that relates unambiguously to some physical earthquake parameter, the seismic moment \(M_0\). Estimating the moment magnitude \(M_w\) requires to observe the low-frequency part of the radiated waveform spectrum, which is often challenging for small events for which low frequencies are under the noise level. In this notebook, we will see how to use tools from the BPMF.spectrum module to estimate moment and local magnitudes.

Moment Magnitude \(M_w\)

For events that were accurately relocated (\(h_{max,unc} < 5\) km) with NonLinLoc, we can try to estimate their seismic moment \(M_0\) and moment magnitude \(M_w\) by fitting their network-averaged, path-corrected displacement spectrum with the omega-square model (Brune or Boatwright model).

We recall that (note that \(\log\) is the base 10 logarithm):

\[M_w = \frac{2}{3} \left( \log M_0 - 9.1 \right)\quad (1),\]

with \(M_0\) in N.m.

Models of displacement spectra in the far-field predict the shape:

\[S(f) = \dfrac{\Omega_0}{\left( 1 + (f/f_c)^{\gamma n} \right)^{1/\gamma}} \exp \left( - \dfrac{\pi \tau^{P/S}(x, \xi) f}{Q^{P/S}} \right)\quad (2).\]

Equation (2) is a general expression where \(f_c\) is the corner frequency and \(n\) controls the high-frequency fall-off rate of the spectrum (typically \(n \approx 2\)) and \(\gamma\) controls the sharpness of the spectrum’s corner. The pre-exponential term is the source term and the exponential term is attenuation. \(\tau^{P/S}(x, \xi)\) is the P/S travel time from source location \(\xi\) to receiver location \(x\), and \(Q^{P/S}\) is the P/S attenuation factor. Two special cases of (2) are: - The Brune spectrum (\(n=2\), \(\gamma = 1\)):

\[S_{\mathrm{Brune}}(f) = \dfrac{\Omega_0}{1 + (f/f_c)^{2}} \exp \left( - \dfrac{\pi \tau^{P/S}(x, \xi) f}{Q^{P/S}} \right)\quad (3)\]

• Boatwright’s spectrum (\(n=2\), \(\gamma=2\)):

\[S_{\mathrm{Boatwright}}(f) = \dfrac{\Omega_0}{\left( 1 + (f/f_c)^{2 n} \right)^{1/2}} \exp \left( - \dfrac{\pi \tau^{P/S}(x, \xi) f}{Q^{P/S}} \right)\quad (4)\]

The low-frequency plateau \(\Omega_0\) is proportional to the seismic moment \(M_0\).

\[\Omega^{P/S}_0 = \dfrac{M_0 R^{P/S}}{4 \pi \sqrt{\rho(\xi) \rho(x) v_{P/S}(\xi) v_{P/S}(x) } v_{P/S}^2(\xi) \mathcal{R}_{P/S}(x, \xi)}\quad (5).\]

In Equation (5), the P/S superscript is for P and S wave, \(x\) is the receiver location and \(\xi\) is the source location. \(R^{P/S}\) is the average radiation pattern (\(R^P = \sqrt{4/15}\) and \(R^S = \sqrt{2/5}\), see Aki and Richards, 2002), \(\rho\) and \(v_{P/S}\) are the medium density and P/S-wave velocity at the source, respectively, \(\mathcal{R}_{P/S}(x, \xi)\) is the P/S ray length (which can be approximated by \(r\), the source-receiver distance). More details in Aki and Richards, 2002 and Boatwright, 1978.

Local Magnitude \(M_L\)

The local magnitude was initially introduced by Richter as the log ratio of peak displacement amplitudes between two earthquakes:

\[M_L = \log \left( \dfrac{A(X)}{A_0(X)} \right),\quad (6)\]

which results in an empirical formula correcting for the source-receiver epicentral distance \(X\):

\[M_L = \log A(X) + \alpha \log X + \beta,\quad (7)\]

where \(\alpha\) and \(\beta\) and region-specific parameters (see Peter Shearer’s Introduction to Seismology book, Equations 9.68 and 9.69).

In this notebook, we define the local magnitude in a slightly different way as:

\[M_L = \log A(r) + \log r + C,\quad (8)\]

where \(r\) is the source-receiver hypocentral distance and \(C\) is a calibration constant. In Equation (8), attenuation is neglected. For small earthquakes and low noise seismograms, \(M_L \propto \log M_0\).

References

Aki, K., & Richards, P. G. (2002). Quantitative seismology.

Boatwright, J. (1978). Detailed spectral analysis of two small New York State earthquakes. Bulletin of the Seismological Society of America, 68(4), 1117-1131.

Shearer, Peter M. Introduction to seismology. Cambridge university press, 2019.

import os
# choose the number of threads you want to limit the computation to
n_CPUs = 24
os.environ["OMP_NUM_THREADS"] = str(n_CPUs)

import BPMF
import h5py as h5
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from BPMF.data_reader_examples import data_reader_mseed
from time import time as give_time

INPUT_DB_FILENAME = "final_catalog.h5"
NETWORK_FILENAME = "network.csv"

net = BPMF.dataset.Network(NETWORK_FILENAME)
net.read()

Load the metadata of the previously detected events

This is similar to what we did in 6_relocate.

events = []
with h5.File(os.path.join(BPMF.cfg.OUTPUT_PATH, INPUT_DB_FILENAME), mode="r") as fin:
    for group_id in fin.keys():
        events.append(
            BPMF.dataset.Event.read_from_file(
                INPUT_DB_FILENAME,
                db_path=BPMF.cfg.OUTPUT_PATH,
                gid=group_id,
                data_reader=data_reader_mseed
            )
        )
        # set the source-receiver distance as this is an important parameter
        # to know to correct for geometrical spreading
        events[-1].set_source_receiver_dist(net)
        if not hasattr(events[-1], "arrival_times"):
            events[-1].set_arrival_times_from_moveouts(verbose=0)
        # set the default moveouts to the theoretical times, that is, those computed in the velocity models
        events[-1].set_moveouts_to_theoretical_times()
        # when available, set the moveouts to the values defined by the empirical (PhaseNet) picks
        events[-1].set_moveouts_to_empirical_times()

print(f"Loaded {len(events)} events.")

Loaded 66 events.

for i, ev in enumerate(events):
    print(i, ev.origin_time)

0 2012-07-26T00:58:16.320000Z
1 2012-07-26T01:15:54.240000Z
2 2012-07-26T01:15:54.320000Z
3 2012-07-26T01:16:29.800000Z
4 2012-07-26T01:16:54.440000Z
5 2012-07-26T01:18:32.720000Z
6 2012-07-26T01:39:51.200000Z
7 2012-07-26T01:47:30.280000Z
8 2012-07-26T02:22:49.640000Z
9 2012-07-26T01:02:53.520000Z
10 2012-07-26T03:00:38.960000Z
11 2012-07-26T03:32:51.000000Z
12 2012-07-26T03:38:13.400000Z
13 2012-07-26T01:10:22.040000Z
14 2012-07-26T01:15:15.880000Z
15 2012-07-26T13:35:26.320000Z
16 2012-07-26T09:21:39.640000Z
17 2012-07-26T09:28:48.560000Z
18 2012-07-26T10:07:23.800000Z
19 2012-07-26T16:27:17.800000Z
20 2012-07-26T11:55:35.080000Z
21 2012-07-26T00:58:11.080000Z
22 2012-07-26T01:15:15.200000Z
23 2012-07-26T01:22:37.880000Z
24 2012-07-26T01:42:22.800000Z
25 2012-07-26T01:52:36.960000Z
26 2012-07-26T00:59:12.800000Z
27 2012-07-26T02:24:36.440000Z
28 2012-07-26T03:05:35.080000Z
29 2012-07-26T04:45:03.960000Z
30 2012-07-26T13:48:32.800000Z
31 2012-07-26T13:50:18.560000Z
32 2012-07-26T13:51:58.560000Z
33 2012-07-26T13:53:31.640000Z
34 2012-07-26T13:54:54.760000Z
35 2012-07-26T13:56:54.240000Z
36 2012-07-26T01:03:47.000000Z
37 2012-07-26T01:09:16.320000Z
38 2012-07-26T01:09:25.960000Z
39 2012-07-26T01:10:44.200000Z
40 2012-07-26T01:12:32.960000Z
41 2012-07-26T15:06:20.040000Z
42 2012-07-26T17:28:21.000000Z
43 2012-07-26T00:20:47.480000Z
44 2012-07-26T01:58:19.800000Z
45 2012-07-26T02:35:01.160000Z
46 2012-07-26T01:35:14.720000Z
47 2012-07-26T14:38:50.760000Z
48 2012-07-26T04:43:38.640000Z
49 2012-07-26T04:46:49.480000Z
50 2012-07-26T04:51:06.800000Z
51 2012-07-26T04:48:38.400000Z
52 2012-07-26T14:49:28.320000Z
53 2012-07-26T03:08:33.640000Z
54 2012-07-26T03:10:54.920000Z
55 2012-07-26T04:30:33.400000Z
56 2012-07-26T05:22:04.320000Z
57 2012-07-26T05:45:45.880000Z
58 2012-07-26T05:46:59.080000Z
59 2012-07-26T05:57:16.320000Z
60 2012-07-26T08:07:38.680000Z
61 2012-07-26T08:08:25.680000Z
62 2012-07-26T10:16:45.640000Z
63 2012-07-26T10:53:46.440000Z
64 2012-07-26T11:02:24.200000Z
65 2012-07-26T01:52:27.800000Z

Example with a single event

In this section, we will demonstrate the workflow we adopt to compute the moment magnitude of each event. It consists in: - Extracting short windows around the P and S waves, as well as a noise window taken before the P wave. - Computing the velocity noise/P/S spectra on each channel. - Correcting for the path effects (geometrical spreading and attenuation). - Averaging the velocity over the network and integrating to get the displacement spectrum. - Fitting the displacement spectrum with a model of our choice (here, the Boatwright model).

#         waveform extraction parameters
# PHASE_ON_COMP: dictionary defining which moveout we use to extract the waveform
PHASE_ON_COMP_S = {"N": "S", "1": "S", "E": "S", "2": "S", "Z": "S"}
PHASE_ON_COMP_P = {"N": "P", "1": "P", "E": "P", "2": "P", "Z": "P"}
# OFFSET_PHASE: dictionary defining the time offset taken before a given phase
#               for example OFFSET_PHASE["P"] = 1.0 means that we extract the window
#               1 second before the predicted P arrival time
OFFSET_PHASE = {"P": 0.5, "S": 0.5}
DURATION_SEC = 3.0
OFFSET_OT_SEC_NOISE = DURATION_SEC + 1.0
TIME_SHIFTED = True
DATA_FOLDER = "raw"
DATA_READER = data_reader_mseed
ATTACH_RESPONSE = True

#       spectral inversion parameters
SPECTRAL_MODEL = "boatwright"
RHO_SOURCE_KGM3 = 2700.0
VS_SOURCE_MS = 3500.0
VP_SOURCE_MS = VS_SOURCE_MS * 1.72
RHO_RECEIVER_KGM3 = 2600.0
VS_RECEIVER_MS = 2800.0
VP_RECEIVER_MS = VS_RECEIVER_MS * 1.72
FREQ_MIN_HZ = 0.5
FREQ_MAX_HZ = 20.
NUM_FREQS = 100
SNR_THRESHOLD = 10.
MIN_NUM_VALID_CHANNELS_PER_FREQ_BIN = 5
MIN_FRACTION_VALID_POINTS_BELOW_FC = 0.20
MAX_RELATIVE_DISTANCE_ERR_PCT = 33.
NUM_CHANNEL_WEIGHTED_FIT = True

EVENT_IDX = 13
event = events[EVENT_IDX]
print(f"The maximum horizontal location uncertainty of event {EVENT_IDX} is {event.hmax_unc:.2f}km.")
print(f"The minimum horizontal location uncertainty of event {EVENT_IDX} is {event.hmin_unc:.2f}km.")
print(f"The maximum vertical location uncertainty is {event.vmax_unc:.2f}km.")

The maximum horizontal location uncertainty of event 13 is 2.12km.
The minimum horizontal location uncertainty of event 13 is 1.73km.
The maximum vertical location uncertainty is 3.00km.

In the following, we define 3 obspy.Stream instances with the noise, P-wave and S-wave traces. Here, it is important to identify the clipped waveforms and set them to zero so that we don’t use them in the future. The identification of clipped waveforms is based on a simple kurtosis test (kurtosis is low for clipped waveforms because the extremum values are often reached).

#                 extract waveforms
# first, read short extract before signal as an estimate of noise
event.read_waveforms(
    DURATION_SEC,
    time_shifted=False,
    data_folder=DATA_FOLDER,
    offset_ot=OFFSET_OT_SEC_NOISE,
    attach_response=ATTACH_RESPONSE,
)
noise = event.traces.copy()
noise.remove_sensitivity()

# then, read signal
event.read_waveforms(
    DURATION_SEC,
    phase_on_comp=PHASE_ON_COMP_P,
    offset_phase=OFFSET_PHASE,
    time_shifted=TIME_SHIFTED,
    data_folder=DATA_FOLDER,
    attach_response=ATTACH_RESPONSE,
)
event.traces.remove_sensitivity()
event.zero_out_clipped_waveforms(kurtosis_threshold=-1)
p_wave = event.traces.copy()

event.read_waveforms(
    DURATION_SEC,
    phase_on_comp=PHASE_ON_COMP_S,
    offset_phase=OFFSET_PHASE,
    time_shifted=TIME_SHIFTED,
    data_folder=DATA_FOLDER,
    attach_response=ATTACH_RESPONSE,
)
event.traces.remove_sensitivity()
event.zero_out_clipped_waveforms(kurtosis_threshold=-1)
s_wave = event.traces.copy()

We can now compute the velocity spectra on each channel.

spectrum = BPMF.spectrum.Spectrum(event=event)
spectrum.compute_spectrum(noise, "noise")
spectrum.compute_spectrum(p_wave, "p")
spectrum.compute_spectrum(s_wave, "s")

spectrum.set_target_frequencies(FREQ_MIN_HZ, FREQ_MAX_HZ, NUM_FREQS)
spectrum.resample(spectrum.frequencies, spectrum.phases)
spectrum.compute_signal_to_noise_ratio("p")
spectrum.compute_signal_to_noise_ratio("s")

Let’s plot the velocity spectra and the signal-to-noise ratio (SNR) spectra.

phase_to_plot = "s"
fig = spectrum.plot_spectrum(phase_to_plot, plot_snr=True)

Before correcting for geometrical spreading, we need to build an attenuation model. In general, body-wave attenuation in the lithosphere is described by he quality factor \(Q\), which is related to frequency by (Aki, 1980):

\[Q \propto f^{n},\ 0.5 \leq n \leq 0.8\quad (6)\]

To build a simple model of attenuation, we assume that \(Q^P = Q^S = Q_0 f^n\). We search for an adequate value of \(Q_0\) and \(n\) in the literature. \(Q_0 = 33\) and \(n = 0.75\) reproduce reasonably well the curves shown in Izgi et al., 2020.

References:

Aki, K. (1980). Attenuation of shear-waves in the lithosphere for frequencies from 0.05 to 25 Hz. Physics of the Earth and Planetary Interiors, 21(1), 50-60.

Izgi, G., Eken, T., Gaebler, P., Eulenfeld, T., & Taymaz, T. (2020). Crustal seismic attenuation parameters in the western region of the North Anatolian Fault Zone. Journal of Geodynamics, 134, 101694.

Q_1Hz = 33.
n = 0.75
Q = Q_1Hz * np.power(spectrum.frequencies, n)

spectrum.set_Q_model(Q, spectrum.frequencies)
spectrum.compute_correction_factor(
    RHO_SOURCE_KGM3, RHO_RECEIVER_KGM3,
    VP_SOURCE_MS, VP_RECEIVER_MS,
    VS_SOURCE_MS, VS_RECEIVER_MS
)

source_parameters = {}
for phase_for_mag in ["p", "s"]:
    spectrum.compute_network_average_spectrum(
        phase_for_mag,
        SNR_THRESHOLD,
        min_num_valid_channels_per_freq_bin=MIN_NUM_VALID_CHANNELS_PER_FREQ_BIN,
        max_relative_distance_err_pct=MAX_RELATIVE_DISTANCE_ERR_PCT,
        verbose=1
        )
    spectrum.integrate(phase_for_mag, average=True)
    spectrum.fit_average_spectrum(
            phase_for_mag,
            model=SPECTRAL_MODEL,
            min_fraction_valid_points_below_fc=MIN_FRACTION_VALID_POINTS_BELOW_FC,
            weighted=NUM_CHANNEL_WEIGHTED_FIT,
            )
    if spectrum.inversion_success:
        rel_M0_err = 100.*spectrum.M0_err/spectrum.M0
        rel_fc_err = 100.*spectrum.fc_err/spectrum.fc
        if rel_M0_err > 10. or spectrum.fc < 0. or spectrum.fc > 25.:
            continue
        print(f"Relative error on M0: {rel_M0_err:.2f}%")
        print(f"Relative error on fc: {rel_fc_err:.2f}%")
        # event.set_aux_data({f"Mw_{phase_for_mag}": spectrum.Mw})
        figtitle = (f"{event.origin_time.strftime('%Y-%m-%dT%H:%M:%S')}: "
                    f"{event.latitude:.3f}""\u00b0"
                    f"N, {event.longitude:.3f}""\u00b0"
                    f"E, {event.depth:.1f}km, "
                    r"$\Delta M_0 / M_0=$"f"{rel_M0_err:.1f}%, "
                    r"$\Delta f_c / f_c=$"f"{rel_fc_err:.1f}%")
        source_parameters[f"M0_{phase_for_mag}"] = spectrum.M0
        source_parameters[f"Mw_{phase_for_mag}"] = spectrum.Mw
        source_parameters[f"fc_{phase_for_mag}"] = spectrum.fc
        source_parameters[f"M0_err_{phase_for_mag}"] = spectrum.M0_err
        source_parameters[f"fc_err_{phase_for_mag}"] = spectrum.fc_err
        fig = spectrum.plot_average_spectrum(
                phase_for_mag,
                plot_fit=True,
                figname=f"{phase_for_mag}_spectrum_{EVENT_IDX}",
                figtitle=figtitle,
                figsize=(8, 8),
                plot_std=True,
                plot_num_valid_channels=True,
                )

Source-receiver distance relative error is too high: 33.42
Source-receiver distance relative error is too high: 33.42
Source-receiver distance relative error is too high: 33.42
Relative error on M0: 3.53%
Relative error on fc: 2.68%
Source-receiver distance relative error is too high: 33.42
Source-receiver distance relative error is too high: 33.42
Source-receiver distance relative error is too high: 33.42
Relative error on M0: 4.23%
Relative error on fc: 2.99%

To get a single final estimate of the moment magnitude \(M_w\), we average the P- and S-wave estimates. Thus,

\[M_w = \frac{1}{2} (M_w^P + M_w^S).\quad (7)\]

Using the definition of the moment magnitude, we can derive the following formula for the magnitude error:

\[d M_w = \frac{1}{3} \left( \frac{dM_0^P}{M_0^P} + \frac{dM_0^S}{M_0^S} \right)\]

Mw_exists = False
norm = 0.
Mw = 0.
Mw_err = 0.
for ph in ["p", "s"]:
    if f"Mw_{ph}" in source_parameters:
        Mw += source_parameters[f"Mw_{ph}"]
        Mw_err += 2./3. * source_parameters[f"M0_err_{ph}"]/source_parameters[f"M0_{ph}"]
        norm += 1
        Mw_exists = True
if Mw_exists:
    Mw /= norm
    Mw_err /= norm
    source_parameters["Mw"] = Mw
    source_parameters["Mw_err"] = Mw_err
else:
    Mw = np.nan
    Mw_err = np.nan

print(f"The P-S averaged moment magnitude is {Mw:.2f} +/- {Mw_err:.2f}")
source_parameters["Mw"] = Mw
source_parameters["Mw_err"] = Mw_err
source_parameters

The P-S averaged moment magnitude is 2.98 +/- 0.03

{'M0_p': 80459275651239.78,
 'Mw_p': 3.203717412293662,
 'fc_p': 7.583974223366631,
 'M0_err_p': 2838966352570.8267,
 'fc_err_p': 0.2032629644982776,
 'M0_s': 16648575393989.844,
 'Mw_s': 2.7475847181211535,
 'fc_s': 6.72861662942346,
 'M0_err_s': 704897217114.1648,
 'fc_err_s': 0.2013588102413158,
 'Mw': 2.9756510652074075,
 'Mw_err': 0.02587476826881796}

Repeat the same processing for every event

#         waveform extraction parameters
# PHASE_ON_COMP: dictionary defining which moveout we use to extract the waveform
PHASE_ON_COMP_S = {"N": "S", "1": "S", "E": "S", "2": "S", "Z": "S"}
PHASE_ON_COMP_P = {"N": "P", "1": "P", "E": "P", "2": "P", "Z": "P"}
# OFFSET_PHASE: dictionary defining the time offset taken before a given phase
#               for example OFFSET_PHASE["P"] = 1.0 means that we extract the window
#               1 second before the predicted P arrival time
OFFSET_PHASE = {"P": 0.5, "S": 0.5}
DURATION_SEC = 3.0
OFFSET_OT_SEC_NOISE = DURATION_SEC + 1.0
TIME_SHIFTED = True
DATA_FOLDER = "raw"
DATA_READER = data_reader_mseed
ATTACH_RESPONSE = True

#       spectral inversion parameters
SPECTRAL_MODEL = "boatwright"
RHO_SOURCE_KGM3 = 2700.0
VS_SOURCE_MS = 3500.0
VP_SOURCE_MS = VS_SOURCE_MS * 1.72
RHO_RECEIVER_KGM3 = 2600.0
VS_RECEIVER_MS = 2800.0
VP_RECEIVER_MS = VS_RECEIVER_MS * 1.72
FREQ_MIN_HZ = 0.5
FREQ_MAX_HZ = 20.
NUM_FREQS = 100
SNR_THRESHOLD = 10.
MIN_NUM_VALID_CHANNELS_PER_FREQ_BIN = 5
MIN_FRACTION_VALID_POINTS_BELOW_FC = 0.20
MAX_RELATIVE_DISTANCE_ERR_PCT = 33.
NUM_CHANNEL_WEIGHTED_FIT = True

#                   attenuation model
Q_1Hz = 33.
n = 0.75
Q = Q_1Hz * np.power(spectrum.frequencies, n)

Let’s run the magnitude computation code for every event. If fitting the network-averaged displacement spectrum is not possible because of poor SNR or location, or if the inverted parameters’ errors are too large then the moment magnitude cannot be estimated and the returned value is nan.

for i, event in enumerate(events):
    print("========================")
    print(f"Processing event {i}")

    #                 extract waveforms
    # first, read short extract before signal as an estimate of noise
    event.read_waveforms(
        DURATION_SEC,
        time_shifted=False,
        data_folder=DATA_FOLDER,
        offset_ot=OFFSET_OT_SEC_NOISE,
        attach_response=ATTACH_RESPONSE,
    )
    noise = event.traces.copy()
    noise.remove_sensitivity()

    # then, read signal
    event.read_waveforms(
        DURATION_SEC,
        phase_on_comp=PHASE_ON_COMP_P,
        offset_phase=OFFSET_PHASE,
        time_shifted=TIME_SHIFTED,
        data_folder=DATA_FOLDER,
        attach_response=ATTACH_RESPONSE,
    )
    event.traces.remove_sensitivity()
    event.zero_out_clipped_waveforms(kurtosis_threshold=-1)
    p_wave = event.traces.copy()

    event.read_waveforms(
        DURATION_SEC,
        phase_on_comp=PHASE_ON_COMP_S,
        offset_phase=OFFSET_PHASE,
        time_shifted=TIME_SHIFTED,
        data_folder=DATA_FOLDER,
        attach_response=ATTACH_RESPONSE,
    )
    event.traces.remove_sensitivity()
    event.zero_out_clipped_waveforms(kurtosis_threshold=-1)
    s_wave = event.traces.copy()

    spectrum = BPMF.spectrum.Spectrum(event=event)
    spectrum.compute_spectrum(noise, "noise")
    spectrum.compute_spectrum(p_wave, "p")
    spectrum.compute_spectrum(s_wave, "s")

    spectrum.set_target_frequencies(FREQ_MIN_HZ, FREQ_MAX_HZ, NUM_FREQS)
    spectrum.resample(spectrum.frequencies, spectrum.phases)
    spectrum.compute_signal_to_noise_ratio("p")
    spectrum.compute_signal_to_noise_ratio("s")

    spectrum.set_Q_model(Q, spectrum.frequencies)
    spectrum.compute_correction_factor(
        RHO_SOURCE_KGM3, RHO_RECEIVER_KGM3,
        VP_SOURCE_MS, VP_RECEIVER_MS,
        VS_SOURCE_MS, VS_RECEIVER_MS
    )

    source_parameters = {}
    for phase_for_mag in ["p", "s"]:
        spectrum.compute_network_average_spectrum(
            phase_for_mag,
            SNR_THRESHOLD,
            min_num_valid_channels_per_freq_bin=MIN_NUM_VALID_CHANNELS_PER_FREQ_BIN,
            max_relative_distance_err_pct=MAX_RELATIVE_DISTANCE_ERR_PCT
            )
        if not phase_for_mag in spectrum.average_spectra:
            continue
        spectrum.integrate(phase_for_mag, average=True)
        spectrum.fit_average_spectrum(
                phase_for_mag,
                model=SPECTRAL_MODEL,
                min_fraction_valid_points_below_fc=MIN_FRACTION_VALID_POINTS_BELOW_FC,
                weighted=NUM_CHANNEL_WEIGHTED_FIT
                )
        if spectrum.inversion_success:
            rel_M0_err = 100.*spectrum.M0_err/spectrum.M0
            rel_fc_err = 100.*spectrum.fc_err/spectrum.fc
            if rel_M0_err > 10. or spectrum.fc < 0. or spectrum.fc > 25.:
                continue
            # event.set_aux_data({f"Mw_{phase_for_mag}": spectrum.Mw})
            figtitle = (f"{event.origin_time.strftime('%Y-%m-%dT%H:%M:%S')}: "
                        f"{event.latitude:.3f}""\u00b0"
                        f"N, {event.longitude:.3f}""\u00b0"
                        f"E, {event.depth:.1f}km, "
                        r"$\Delta M_0 / M_0=$"f"{rel_M0_err:.1f}%, "
                        r"$\Delta f_c / f_c=$"f"{rel_fc_err:.1f}%")
            source_parameters[f"M0_{phase_for_mag}"] = spectrum.M0
            source_parameters[f"Mw_{phase_for_mag}"] = spectrum.Mw
            source_parameters[f"fc_{phase_for_mag}"] = spectrum.fc
            source_parameters[f"M0_err_{phase_for_mag}"] = spectrum.M0_err
            source_parameters[f"fc_err_{phase_for_mag}"] = spectrum.fc_err
            fig = spectrum.plot_average_spectrum(
                    phase_for_mag,
                    plot_fit=True,
                    figname=f"{phase_for_mag}_spectrum_{i}",
                    figtitle=figtitle,
                    figsize=(8, 8),
                    plot_std=True,
                    plot_num_valid_channels=True,
                    )

    Mw_exists = False
    norm = 0.
    Mw = 0.
    Mw_err = 0.
    for ph in ["p", "s"]:
        if f"Mw_{ph}" in source_parameters:
            Mw += source_parameters[f"Mw_{ph}"]
            Mw_err += 2./3. * source_parameters[f"M0_err_{ph}"]/source_parameters[f"M0_{ph}"]
            norm += 1
            Mw_exists = True
    if Mw_exists:
        Mw /= norm
        Mw_err /= norm
        source_parameters["Mw"] = Mw
        source_parameters["Mw_err"] = Mw_err
    else:
        Mw = np.nan
        Mw_err = np.nan

    print(f"The P-S averaged moment magnitude is {Mw:.2f} +/- {Mw_err:.2f}")
    source_parameters["Mw"] = Mw
    source_parameters["Mw_err"] = Mw_err

    # save all this new information in BPMF.dataset.Event.aux_data
    event.set_aux_data(source_parameters)

========================
Processing event 0
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 1
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 2
Spectrum is below SNR threshold everywhere, cannot fit it.
Spectrum is below SNR threshold everywhere, cannot fit it.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 3
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 4
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 5
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 6
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 7
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 8
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 9
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 10
Spectrum is below SNR threshold everywhere, cannot fit it.
Spectrum is below SNR threshold everywhere, cannot fit it.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 11
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 12
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 13
The P-S averaged moment magnitude is 2.98 +/- 0.03
========================
Processing event 14
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 15
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 16
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 17
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 18
Spectrum is below SNR threshold everywhere, cannot fit it.
Spectrum is below SNR threshold everywhere, cannot fit it.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 19
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 20
Spectrum is below SNR threshold everywhere, cannot fit it.
Spectrum is below SNR threshold everywhere, cannot fit it.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 21
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 22
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 23
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 24
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 25
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 26
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 27
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 28
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 29
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 30
The P-S averaged moment magnitude is 2.92 +/- 0.03
========================
Processing event 31
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 32
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 33
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 34
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 35
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 36
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 37
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 38
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 39
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 40
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 41
Spectrum is below SNR threshold everywhere, cannot fit it.
Spectrum is below SNR threshold everywhere, cannot fit it.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 42
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 43
Spectrum is below SNR threshold everywhere, cannot fit it.
Spectrum is below SNR threshold everywhere, cannot fit it.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 44
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 45
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 46
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 47
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 48
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 49
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 50
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 51
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 52
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 53
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 54
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 55
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 56
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 57
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 58
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 59
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 60
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 61
Spectrum is below SNR threshold everywhere, cannot fit it.
Spectrum is below SNR threshold everywhere, cannot fit it.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 62
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 63
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 64
The P-S averaged moment magnitude is nan +/- nan
========================
Processing event 65
Class instance does not have a `cov_mat` attribute.
The P-S averaged moment magnitude is nan +/- nan

# for event in backprojection_events:
#     event.update_aux_data_database()
BPMF.utils.donefun()

⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
    ⠀⠀⠀⠀⠀⠀⢀⡤⣤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡀⠀⠀⠀⠀⠀⠀
    ⠀⠀⠀⠀⠀⢀⡏⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⣀⠴⠋⠉⠉⡆⠀⠀⠀⠀⠀
    ⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠈⠉⠉⠙⠓⠚⠁⠀⠀⠀⠀⣿⠀⠀⠀⠀⠀
    ⠀⠀⠀⠀⢀⠞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠹⣄⠀⠀⠀⠀
    ⠀⠀⠀⠀⡞⠀⠀⠀⠀⠀⠶⠀⠀⠀⠀⠀⠀⠦⠀⠀⠀⠀⠀⠸⡆⠀⠀⠀
    ⢠⣤⣶⣾⣧⣤⣤⣀⡀⠀⠀⠀⠀⠈⠀⠀⠀⢀⡤⠴⠶⠤⢤⡀⣧⣀⣀⠀
    ⠻⠶⣾⠁⠀⠀⠀⠀⠙⣆⠀⠀⠀⠀⠀⠀⣰⠋⠀⠀⠀⠀⠀⢹⣿⣭⣽⠇
    ⠀⠀⠙⠤⠴⢤⡤⠤⠤⠋⠉⠉⠉⠉⠉⠉⠉⠳⠖⠦⠤⠶⠦⠞⠁⠀⠀⠀
                ALL DONE!⠀⠀⠀⠀

To be continued…

Work-in-progress!