Enriching the GEOFON seismic catalog with automatic energy magnitude estimations

Abstract

We present a seismic catalog (, https://doi.org/10.5880/GFZ.2.6.2023.010) including energy magnitude $M_{\mathrm{e}}$ estimated from P waves recorded at teleseismic distances in the range $\mathrm{20}\mathit{°}\le \mathrm{\Delta }\le \mathrm{98}\mathit{°}$ and for depths shorter than 80 km. The catalog is built starting from the event catalog disseminated by GEOFON (GEOFOrschungsNetz), considering 6349 earthquakes with moment magnitude $M_{\mathrm{w}}\ge \mathrm{5}$ occurring between 2011 and 2023. Magnitudes are computed using 1 031 396 freely available waveforms archived in EIDA (European Integrated Data Archive) and IRIS (Incorporated Research Institutions for Seismology) repositories, retrieved through the standard International Federation of Digital Seismograph Networks (FDSN) web services (https://www.fdsn.org/webservices/, last access: March 2024). A reduced, high-quality catalog for events with $M_{\mathrm{w}}\ge \mathrm{5.8}$ and from which stations and events with only few recordings were removed forms the basis of a detailed analysis of the residuals of individual station measurements, which are decomposed into station- and event-specific terms and a term accounting for remaining variability. The derived $M_{\mathrm{e}}$ values are compared to $M_{\mathrm{w}}$ computed by GEOFON and with the $M_{\mathrm{e}}$ values calculated by IRIS. Software and tools developed for downloading and processing waveforms for bulk analysis and an add-on for SeisComP for real-time assessment of $M_{\mathrm{e}}$ in a monitoring context are also provided alongside the catalog. The SeisComP add-on has been part of the GEOFON routine processing since December 2021 to compute and disseminate $M_{\mathrm{e}}$ for major events via the existing services.

1 Introduction

Several magnitude scales have been defined to characterize the size of an earthquake. We can, however, divide magnitude scales in two groups, with one including magnitudes based on the amplitudes and periods of different seismic phases measured on band-limited signals (e.g., the body and surface wave magnitudes; ) and the other including magnitude scales related to estimations of macroscopic physical parameters of the earthquake source. The latter comprise the moment ($M_{\mathrm{w}}$; ) and the energy ($M_{\mathrm{e}}$; ) magnitudes, which are based on seismic moment and radiated seismic energy , respectively. These two magnitude scales are somewhat complementary because, although both represent an estimation of earthquake-related energy, they are determined by different parts of the source spectrum. The seismic moment extrapolated from the low-frequency end and represents the release of elastic energy stored in the Earth's crust or mantle being proportional to the integrated slip across the fault surface. The radiated seismic energy describes the fraction of the total energy released being radiated as seismic waves across all frequencies; i.e., it depends on the earthquake dynamics such as rupture velocity but also stress drop. $M_{\mathrm{e}}$ estimates have been shown to play an important role when used in conjunction with $M_{\mathrm{w}}$ to better characterize the tsunami and shaking potential of an earthquake .

$M_{\mathrm{w}}$ is routinely computed from long period signals of broadband recordings, and it has become a robust and reliable source parameter for large and moderate earthquakes worldwide . On the other hand, the computation of $M_{\mathrm{e}}$ is hindered by the necessity of integrating the velocity power spectra over a wide frequency range, while using signals in a limited bandwidth and taking into account propagation effects at high frequencies.

Aiming at validating and testing the procedures for operational purposes, we present a seismic catalog of $M_{\mathrm{e}}$ computed, following the methodology proposed by and for the rapid assessment of energy magnitude (i.e., without requiring additional source information other than the hypocentral location). The approach is based on the analysis of spectra computed for teleseismic vertical component P waveforms. Teleseismic P waves are commonly used to compute $M_{\mathrm{e}}$ for global earthquakes, as their energy loss during propagation can be more reliably modeled compared to S waves. We further present a detailed analysis of the residuals in a reduced high-quality catalog for events with $M_{\mathrm{w}}\ge \mathrm{5.8}$ with respect to the $M_{\mathrm{w}}$ available in the GEOFON (GEOFOrschungsNetz) catalog and the $M_{\mathrm{e}}$ values computed by IRIS.

2 Energy magnitude computation

2.1 Single-station estimation

We implement the methodology proposed by and to compute $M_{\mathrm{e}}$. Teleseismic vertical component P waveforms (BHZ channels, where B stands for broadband, H for high-gain seismometer, and Z for vertical component) are analyzed in the distance range from $\mathrm{20}\mathit{°}$ to $\mathrm{98}\mathit{°}$ and for earthquakes shallower than 80 km. Standard teleseismic range usually starts at $\mathrm{30}\mathit{°}$, but we use $\mathrm{20}\mathit{°}$ to allow closer stations to be used for rapid response purposes. Shorter distances, however, are difficult to include for global earthquakes, as regional effects are not well accounted for with a 1-D model. Propagation effects are accounted for by frequency-dependent amplitude decay functions that are computed numerically for the ak135Q model in the frequency range 0.012–1 Hz.

An estimate of radiated seismic energy $E_{\mathrm{s}}$ is obtained for each single station from the integral of the power spectra of the vertical component P waveform, which is corrected for propagation effects :

1$$E_{\mathrm{s}}=\left[{\displaystyle \frac{\mathrm{2}}{\mathrm{15}\mathit{\pi }\mathit{\rho }{\mathit{\alpha }}^{\mathrm{5}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{5}\mathit{\pi }\mathit{\rho }{\mathit{\beta }}^{\mathrm{5}}}}\right]\underset{f_{\mathrm{1}}}{\overset{f_{\mathrm{2}}}{\int }}{\left|{\displaystyle \frac{\dot{u}\left(f\right)}{G\left(f\right)/\mathrm{2}\mathit{\pi }f}}\right|}^{\mathrm{2}}\mathrm{d}f,$$ where $\mathit{\alpha }$, $\mathit{\beta }$, and $\mathit{\rho }$ are the P-wave velocity, S-wave velocity, and the density at the source, respectively. $f$ is the frequency, and $f_{\mathrm{1}}=\mathrm{0.012}$ Hz and $f_{\mathrm{2}}=\mathrm{1}$ Hz are the lower and upper limits of the considered spectral bandwidth. $\dot{u}\left(f\right)$ is the P-wave velocity spectrum. $G\left(f\right)$ is the median value of Green's function spectrum for displacement, computed from multiple combinations of focal mechanisms and varying strike, dip, and rake over a regular grid .

We used analysis windows starting 10 s before the P arrival and with lengths of 90 s for $M_{\mathrm{w}}\le \mathrm{7.5}$, 120 s for $\mathrm{7.5}<M_{\mathrm{w}}\le \mathrm{8.5}$, and 180 s for $M_{\mathrm{w}}>\mathrm{8.5}$. The energy magnitude $M_{\mathrm{e}}$ estimate for a single-event station pair is in turn computed as $M_{\mathrm{e}}=\mathrm{2}/\mathrm{3}({\log }_{\mathrm{10}}{E}_{\mathrm{s}}-\mathrm{4.4})$, with $E_{\mathrm{s}}$ given in joules . The procedure provides $M_{\mathrm{e}}$ estimates at each recording station that can be averaged to minimize path-specific deviations not accounted for by the theoretical model (e.g., directivity and focal mechanism effects and regional variations in attenuation).

2.2

Open-source tool for computing $M_{\mathrm{e}}$

The above procedure is implemented in the package me-compute . The program uses stream2segment to download events, station metadata, and waveforms from International Federation of Digital Seismograph Networks (FDSN)-compliant repositories in a structured query language (SQL) database.

In our application, the download is configured to fetch events from the GEOFON event web service, selecting events with computed $M_{\mathrm{w}}$ in the time span 2011–2023. Waveforms are download from EIDA and IRIS (https://service.iris.edu/, last access: March 2024) data centers. The processing routine is implemented in a Python module, which computes the station energy magnitude for each downloaded waveform segment, as summarized in Sect. , and then calculates the event energy magnitude $M_{\mathrm{e}}$ as the mean of all station magnitudes within the 5th–95th percentile range.

The final output consists of the following files:

• a tabular file in a hierarchical data format (HDF), where each row represents the metadata and measurements, specifically also the station energy magnitude estimate, for a single waveform;

• a tabular file in a comma-separated value (CSV) format aggregating the results of the previous file, where each row represents a seismic event, reporting the event data end metadata and including the $M_{\mathrm{e}}$ estimate for the event;

• an HTML (HyperText Markup Language) file visualizing the selected content reported in the CSV file, where the information for each event can be visualized on an interactive map; and

• one file per processed event in the Quake Markup Language (QuakeML) format, where we also included the $M_{\mathrm{e}}$ value.

All files produced by me-compute are disseminated in the data archive (; 10.5880/GFZ.2.6.2023.010), along with the stream2segment and me-compute configuration files.

3 Catalog compilation

We use me-compute to compute $M_{\mathrm{e}}$ for $M_{\mathrm{w}}\ge \mathrm{5}$ earthquakes since 2011 in the GEOFON catalog. Table 1 summarizes the steps followed to compile the disseminated $M_{\mathrm{e}}$ catalog. The catalog reports the single-waveform energy magnitude $M_{\mathrm{e}ij}$ estimated at station $j$ for earthquake $i$. The energy magnitude $M_{\mathrm{e}}$ for each considered event $i$ is then computed as the median of $M_{\mathrm{e}ij}$ over the set of recording stations, without considering station static corrections. The starting dataset D0 consists of more than 1 000 000 waveforms (channels BHZ) generated by 6963 earthquakes recorded by 7765 stations belonging to 246 different networks. Only recordings with an average signal-to-noise ratio (SNR) for the amplitude greater than 3 within the frequency range of interest are included in D0. Several integrity and quality checks are applied to remove outliers and faulty signals. Dataset D1 is obtained by analyzing the median residual at the network level, discarding 14 networks characterized by median residuals outside the 2.5 and 97.5 percentile range (Fig. a). Dataset D2 is then generated by analyzing the station median values and excluding 382 stations with residuals outside the 2.5–97.5 percentile range (Fig. b). Most of the networks and stations removed will have instrumental problems or faulty metadata regarding instrument responses, although in some cases stations with very strong site effects might also be excluded.

Figure 1

Median network residuals (circles) for dataset D0 (a) and median station residuals for dataset D1 (b). Red lines correspond to the 2.5 and 97.5 percentiles of the distributions. For each network (a) and station (b), the horizontal bars correspond to the interval median (circle) $\pm \mathrm{1}$ median absolute deviation (MAD). A few values falling outside the range considered for the horizontal axis are not shown.

[Figure omitted. See PDF]

An anomaly score is computed to further refine the dataset by flagging anomalous amplitudes using the software sdaas . The software, developed from the work of , is based on a machine learning algorithm specifically designed for outlier detection (Isolation Forest) which computes an anomaly score in $[\mathrm{0},\mathrm{1}]$, representing the degree of belief of a waveform to be an outlier. The score can be used to assign robustness weights or to define thresholds above which data can be discarded. After inspecting the distribution of the anomaly scores, we set the threshold to 0.62 for $M_{\mathrm{w}}<\mathrm{7.5}$ and to 0.80 for $M_{\mathrm{w}}\ge \mathrm{7.5}$.

The spatial distribution of events and stations generating dataset D3 are shown in Fig. a, b. The corresponding $M_{\mathrm{e}}$ residuals are shown in Fig.  against distance and $M_{\mathrm{w}}$. The largest positive residuals correspond mostly to earthquakes with $M_{\mathrm{w}}<\mathrm{6}$ recorded at distances $\mathrm{\Delta }>\mathrm{60}\mathit{°}$, where the implemented methodology is expected to generate biased station $M_{\mathrm{e}}$ estimates due to the limitations in the analyzed bandwidth and low signal-to-noise ratio . The overall residual distribution is unbiased and does not show trends of the mean value with distance and magnitude.

Figure 2

Panels (a) and (b) show event and station locations for dataset D3 (Table 1), respectively. Panels (c) and (d) show event and station locations for dataset D6 (Table 1), respectively.

[Figure omitted. See PDF]

Table 1

Datasets considered in this study. Datasets in bold are available in .

Dataset	Records	Networks	Stations	Events	Selections applied sequentially
D0	1 126 465	246	7765	6963	$M_{\mathrm{w}}\ge \mathrm{5}$
D1	1 072 381	232	7617	6944	Network selection (2.5–97.5 %)
D2	1 034 833	228	7235	6880	Station selection (2.5–97.5 %)
D3	1 031 396	228	7234	6349	Anomaly score ($<\mathrm{0.62}$, $<\mathrm{0.8}$ for $M_{\mathrm{w}}<\mathrm{7.5},\ge \mathrm{7.5}$, resp.)
D4	754 025	228	7228	1731	$M_{\mathrm{w}}\ge \mathrm{5.8}$
D5	751 567	227	7135	1731	No. of records per station $\ge \mathrm{10}$
D6	750 903	227	7135	1671	No. of records per event $\ge \mathrm{10}$
Dg				153	Comparison between D6 and real time

Figure 3

Energy magnitude residuals versus distance (a) and moment magnitude (b) for dataset D3. Blue dots indicate residuals also included in D6. The horizontal red lines bound the 90 % confidence interval $[-\mathrm{0.43},\mathrm{0.50}]$ of the residual distribution. The error bars indicate the mean $\pm \mathrm{1}$ standard deviation of the residuals computed over different distance ($\mathrm{20}\mathit{°}$ wide) and magnitude (1 m.u. wide) intervals.

[Figure omitted. See PDF]

Therefore, we further limit the dataset by only considering events with $M_{\mathrm{w}}\ge \mathrm{5.8}$ and at least 10 single-station measurements; we further exclude stations with fewer than 10 recordings in total. We added a column in the disseminated D3 dataset to flag lines corresponding to D6. It consists of $\sim \mathrm{750}\mathrm{000}$ waveforms for 1671 earthquakes and 7135 stations. The event and station locations of D6 are shown in Fig. c and d.

4 Quality assessment via residual analysis

We perform residual analysis to validate the D6 catalog. The relationship between $M_{\mathrm{e}}$ and $M_{\mathrm{w}}$ is analyzed by performing the following mixed-effects regression :

2$$M_{\mathrm{e}ij}=c_{\mathrm{1}}+c_{\mathrm{2}}{M}_{\mathrm{w}i}+\mathit{\delta }{S}_j+\mathit{\delta }{E}_i+{\mathit{ϵ}}_{ij},$$ where $M_{\mathrm{e}ij}$ is the single-waveform energy magnitude estimate at station $j$ for earthquake $i$, $M_{\mathrm{w}i}$ is the moment magnitude of earthquake $i$, intercept $c_{\mathrm{1}}$ and slope $c_{\mathrm{2}}$ parameters define the median model, and $\mathit{\delta }{S}_i$ and $\mathit{\delta }{E}_j$ are terms that capture station-specific and earthquake-specific adjustments, respectively. ${\mathit{ϵ}}_{ij}$ accounts for the leftover effects (i.e., residuals that are specific to a particular path/waveform). The random effects $\mathit{\delta }S$, $\mathit{\delta }E$, and $\mathit{ϵ}$ are zero mean normal distributions by construction. In particular, $\mathit{\delta }{S}_j$ (inter-station residual) can represent site effects or instrumental gain corrections, with most of the latter probably removed by the outlier filtering stages described above. The interevent residual $\mathit{\delta }{E}_i$ is an event-specific deviation from the $M_{\mathrm{e}}$ expected for a given $M_{\mathrm{w}}$ from the linear regression term. Finally, ${\mathit{ϵ}}_{ij}$ can be thought of as a noise term for individual measurements, which can be either related to path-specific heterogeneity in attenuation with respect to the 1-D reference model or the influence of ambient noise on the actual measurement.

The interevent and inter-station term distributions are shown in Fig. , which are described by standard deviations of $\mathit{\tau }=\mathrm{0.246}$ and ${\mathit{\varphi }}_{\mathrm{S}}=\mathrm{0.188}$ m.u., respectively; the standard deviation of the $\mathit{ϵ}$ is ${\mathit{\varphi }}_{\mathrm{0}}=\mathrm{0.232}$ m.u. When combining the interevent variability $\mathit{\tau }$ with the intra-event variability equal to $\mathit{\varphi }=\sqrt{{\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{2}}+{\mathit{\varphi }}_{\mathrm{S}}^{\mathrm{2}}}$, we obtain the total standard deviation $\mathit{\sigma }=\sqrt{{\mathit{\tau }}^{\mathrm{2}}+{\mathit{\varphi }}^{\mathrm{2}}}=\mathrm{0.407}$, which represents the variability in the single-station $M_{\mathrm{e}ij}$ residuals with respect to the average $M_{\mathrm{e}}$ computed per event. It is worth noting that the $\mathit{\delta }{S}_j$ values can be used as station corrections to compute the energy magnitude of new events. In this case, the inter-station contribution to the total variability is removed, and the expected variability in the $M_{\mathrm{e}ij}$ distribution is reduced to $\sqrt{{\mathit{\tau }}^{\mathrm{2}}+{\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{2}}}=\mathrm{0.338}$. Finally, the linear regression model is defined by the coefficients $c_{\mathrm{1}}=(\mathrm{0.77}\pm \mathrm{0.09})$ m.u. and $c_{\mathrm{2}}=(\mathrm{0.92}\pm \mathrm{0.01})$. Considering the simplicity of the linear model in Eq. () and the large dataset analyzed, the uncertainty in the median model (sometimes referred to as ${\mathit{\sigma }}_{\mathit{\mu }}$; ) is very low, increasing from 0.007 for $M_{\mathrm{w}}=\mathrm{6}$ to 0.039 for $M_{\mathrm{w}}=\mathrm{9}$.

Figure 4

Cumulative distribution functions for event $\mathit{\delta }E$ (a), station $\mathit{\delta }S$ (b), and leftover $\mathit{ϵ}$ distributions (circles) determined according to the mixed-effects regression in Eq. () applied to dataset D6. Dotted lines correspond to standard deviations $\pm \mathrm{1}\mathit{\tau }$ (a), $\pm \mathrm{1}{\mathit{\varphi }}_S$ (b), and $\pm \mathrm{1}{\mathit{\varphi }}_{\mathrm{0}}$ (c). Horizontal red lines in panels (a) and (b) are the standard errors in the random effects. In panel (c), values of $\mathit{ϵ}$ exceeding $\pm \mathrm{1.2}$ in absolute value are not shown.

[Figure omitted. See PDF]

Figure 5

(a) Distribution of the site-specific residuals $\mathit{\delta }S$ (see Eq. ) and (b) magnified view over a portion of Europe. Numbers in panel (b) indicate the following locations: (1) Netherlands; (2) Harz highlands, Germany; (3) Switzerland; (4) Po plain, Italy; (5) Pyrenees mountain range; (6) Apennine Mountains; (7) East Anatolian Fault region; and (8) Moesian Platform.

[Figure omitted. See PDF]

We show the spatial distribution of $\mathit{\delta }S$ in Fig. . Since $M_{\mathrm{e}ij}$ is computed considering spectral values below 1 Hz and using teleseismic recordings for distances above $\mathrm{20}\mathit{°}$, $\mathit{\delta }S$ capture station-specific effects connected to large-scale geological and tectonic crustal features, as exemplified in Fig. b for stations located in Europe. Positive $\mathit{\delta }S$ (i.e., $M_{\mathrm{e}ij}$ larger than the median) are observed for stations located in basins like in the Po plain, in the Moesian region, in the Netherlands, and in the East Anatolian Fault region. Negative values $\mathit{\delta }S$ (i.e., $M_{\mathrm{e}ij}$ lower than the median) are observed for stations located in mountain ranges such as the Pyrenees, the Alps, or in Harz Highlands and also tectonically highly active regions that are as cratonic as the East African Rift. The station terms can represent both site amplification, e.g., for stations in sedimentary basins, and anomalously high or low attenuation in the crust and or mantle surrounding the station. The station-specific residuals are disseminated along with the catalog to allow the computation of $M_{\mathrm{e}}$ for future earthquakes, while taking into account static magnitude corrections to reduce variability.

The spatial distribution of the interevent variability, $\mathit{\delta }E$, is shown in Fig.  for the smallest and largest values.

Figure 6

Extreme values for event specific residuals $\mathit{\delta }E$ for the $M_{\mathrm{e}ij}$ versus $M_{\mathrm{w}}$ mixed-effects model of Eq. (). Only values below the 10th percentile (panels b and d) and above the 90th percentile (panels a and c) of the distribution are shown (the percentiles are about $\pm \mathrm{0.3}$). In panels (a) and (b), earthquakes with hypocentral depths shallower than 30 km are selected. In panels (c) and (d), events deeper than 30 km are considered. The distribution of $\mathit{\delta }E$ versus depth for all events is shown in panel (e).

[Figure omitted. See PDF]

Considering depths shallower than 30 km (Fig. a and b), continental Asia, the Philippines, and Indonesia, as well as the Aleutian Islands, show positive values. California, Mexico, central America, and the Mid-Atlantic Ridge are characterized mostly by negative values. Considering deeper events (Fig. c and d), Japan and the Philippines have mostly positive values, while Mexico and central America have mostly negative values. The event-specific residuals are also disseminated, along with the catalog to increase the usefulness of the product from the event point of view, to allow the user to perform further refinements.

Figure 7

(a) Residual distribution $\mathit{ϵ}$ of Eq. () for three different receiving areas, showing only $\left|\mathit{ϵ}\right|>\mathrm{0.30}$. (b) As in panel (a) but considering only the European receiving area. (c) As in panel (a) but considering only the receiving area in California. (d) As in panel (a) but considering only the receiving area in Australia. Circles indicate the earthquake locations.

[Figure omitted. See PDF]

Figure 8

$M_{\mathrm{e}ij}$ versus $M_{\mathrm{w}}$ scaling. Gray circles are the station $M_{\mathrm{e}ij}$ estimates, and filled circles represent event $M_{\mathrm{e}}$ values calculated as medians of all station estimates for that event. Color indicates how many stations contributed to each estimate. The best-fit line in green is derived from the mixed-effects regression (Eq. ), considering $\pm \mathrm{1}$ interevent standard deviation $\mathit{\tau }$ (red lines). The faint black line shows equality for reference.

[Figure omitted. See PDF]

Path-specific residuals $\mathit{ϵ}$ are shown in Fig.  for three selected receiving areas in Europe, California, and Australia. Since in the partition of the residuals the leftover distribution $\mathit{ϵ}$ represents the component not related to systematic station and event effects, they are mostly connected to lateral variability in attenuation in the Earth's interior with respect to the used global 1-D model and amplitude variation related to P-wave radiation patterns for different focal mechanisms.

Finally, the $M_{\mathrm{e}ij}$ versus $M_{\mathrm{w}}$ scaling defined by the linear regression coefficients $c_{\mathrm{1}}$ and $c_{\mathrm{2}}$ of Eq. () is shown in Fig. .

4.1 Catalog validation: comparison with IRIS

The energy magnitude computed in this study is compared to the values disseminated by IRIS through the SPUD service . The methodology implemented by IRIS is described by and based on the analysis of and . Similar to our approach, the energy flux is computed from the P-wave group (P $+$ pP $+$ sP) in the frequency domain. The single-station estimations are corrected for frequency-dependent anelastic attenuation effects and converted back to the energy radiated by the source by applying corrections for geometrical spreading, depth, and mechanism-dependent effects for P waves and considering a theoretical partition of the energy between P and S waves. The energy is computed considering the frequency range 0.014–2 Hz (broadband for $M_{\mathrm{e}}\left(\mathrm{BB}\right)$) or 0.5-02 Hz (high frequency for $M_{\mathrm{e}}\left(\mathrm{HF}\right)$), analyzing stations in the distance range 25–80°. The duration of the time window used for the computation is based on analysis of the cumulative high-frequency energy (0.5–2 Hz) as a function of time. The crossover time used to compute the energy flux is identified at the intersection between the near-constant increasing rate for short times and the relative flat asymptotic behavior for long durations. The SPUD service disseminates both the high-frequency $M_{\mathrm{e}}\left(\mathrm{HF}\right)$ and broadband $M_{\mathrm{e}}\left(\mathrm{BB}\right)$ estimates.

Two regression models are calibrated against the broadband and high-frequency estimates disseminated by IRIS through SPUD. The best-fit models, shown in Fig. , are $M_{\mathrm{e}}=(-\mathrm{0.076}\pm \mathrm{0.229})+(\mathrm{1.002}\pm \mathrm{0.033})M_{\mathrm{e}}\left(\mathrm{HF}\right)$ and $M_{\mathrm{e}}=(\mathrm{0.795}\pm \mathrm{0.188})+(\mathrm{0.896}\pm \mathrm{0.027})M_{\mathrm{e}}\left(\mathrm{BB}\right)$, with the standard deviation of the residuals equal to 0.234 and 0.175, respectively. For the magnitude range from 6 to 8, this results in biases of 0.06 m.u. for $M_{\mathrm{e}}$ vs. $M_{\mathrm{e}}\left(\mathrm{HF}\right)$ and varying from 0.17 to $-\mathrm{0.04}$ m.u. for $M_{\mathrm{e}}$ vs. $M_{\mathrm{e}}\left(\mathrm{BB}\right)$; i.e., our estimates are nearly unbiased relative to $M_{\mathrm{e}}\left(\mathrm{HF}\right)$ and tend to slightly overestimate $M_{\mathrm{e}}\left(\mathrm{BB}\right)$ at the lower end of the applicability range.

Figure 9

Comparison with energy magnitude disseminated by IRIS considering (a) $M_{\mathrm{e}}\left(\mathrm{HF}\right)$ and (b) $M_{\mathrm{e}}\left(\mathrm{BB}\right)$ (717 common events). The red line shows the linear regression fit, and the dotted lines show 1 standard deviation of the $M_{\mathrm{e}}$ residuals. The blue line shows line of equality for reference.

[Figure omitted. See PDF]

4.2 Catalog validation: role of style of faulting

The faulting style is classified into normal, reverse, and strike slip categories, based on the plunge of the P, T, and N axes , as extracted from the GEOFON moment tensor solutions. There is a normal fault (NF) if the plunge ($P$) $\ge \mathrm{60}\mathit{°}$, a strike slip (SS) if the plunge ($N$) $\ge \mathrm{60}\mathit{°}$, and a thrust fault (TF) if the plunge $T\ge \mathrm{50}\mathit{°}$. In the other cases, the earthquake is labeled with OF (other faulting styles). To investigate the role of the style of faulting (SOF), we separate the event term into a fixed offset for each SOF class and a perturbation term for each event. If we indicate the classes of the SOF grouping factor (corresponding to NF, SS, TF, and OF) with $k=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}$ and the class of event $i$ with $k_i$, then the equation for the extended mixed-effects model is

3$$M_{\mathrm{e}ij}=e_{\mathrm{1}}+e_{\mathrm{2}}{M}_{\mathrm{w}i}+\mathit{\delta }{S}_j+[\mathit{\delta }{\mathrm{SOF}}_{k_i}+\mathit{\delta }{E}_{{\mathrm{SOF}}_i}]+{\mathit{ϵ}}_{ij},$$ where $\mathit{\delta }\mathrm{SOF}$ are the terms characterizing the average effects of the different SOF and $\mathit{\delta }{E}_{\mathrm{SOF}}$ are accounting for interevent differences within each SOF class (nested random effects). The standard deviations of the $\mathit{\delta }S$, $\mathit{\delta }\mathrm{SOF}$, $\mathit{\delta }{E}_{\mathrm{SOF}}$, and $\mathit{ϵ}$ distributions are ${\mathit{\varphi }}_{\mathrm{S}}=\mathrm{0.190}$, ${\mathit{\tau }}_{\mathrm{SOF}}=\mathrm{0.095}$, $\mathit{\tau }=\mathrm{0.236}$, and ${\mathit{\varphi }}_{\mathrm{0}}=\mathrm{0.232}$, respectively, generating a total standard deviation $\mathit{\sigma }=\mathrm{0.393}$. The SOF terms are $\mathit{\delta }{\mathrm{SOF}}_{\mathrm{1}}=\mathrm{0.098}$ (NF), $\mathit{\delta }{\mathrm{SOF}}_{\mathrm{2}}=-\mathrm{0.108}$ (SS), $\mathit{\delta }{\mathrm{SOF}}_{\mathrm{3}}=-\mathrm{0.045}$ (TF), and $\mathit{\delta }{\mathrm{SOF}}_{\mathrm{4}}=\mathrm{0.055}$ (OF) (Fig. ). The largest difference is between SS and NF, a total of 0.206 m.u. There is a systematic impact of the SOF on the intercept of the model, but the associated variability is smaller compared to the interevent variability $\mathit{\tau }$ (in other words, SOF effects are statistically significant, but the distributions of interevent terms separated according to faulting style are strongly overlapping).

Figure 10

$M_{\mathrm{e}}$ versus $M_{\mathrm{w}}$ categorized with SOF.

[Figure omitted. See PDF]

The SOF effects might arise due to physical differences (on average) between the different faulting types, e.g., due to systematically different stress drops, differences in the maturity of faults or typical environments (intra-plate vs. interplate), where different faulting types occur most often, or they might be artifacts due to the fact that the method used here does not account for radiation pattern effects, and the teleseismic arrivals utilized here preferentially sample certain parts of the focal sphere. Therefore, we also investigate the role of the SOF in the relationship between $M_{\mathrm{e}}$ derived in this study and the $M_{\mathrm{e}}\left(\mathrm{HF}\right)$ and $M_{\mathrm{e}}\left(\mathrm{BB}\right)$ values disseminated by IRIS. We recall that the methodology implemented by IRIS accounts for radiation pattern effects, which are related to the SOF. For this analysis, the regression model is the following: 4$$M_{\mathrm{e}}=g_{\mathrm{1}}+g_{\mathrm{2}}{M}_{\mathrm{iris}}+\mathit{\delta }\mathrm{SOF}+\mathit{ϵ},$$ where $M_{\mathrm{iris}}$ is either $M_{\mathrm{e}}\left(\mathrm{HF}\right)$ or $M_{\mathrm{e}}\left(\mathrm{BB}\right)$. Results shown in Fig.  confirm that the largest intercept difference is between normal and strike slip events, and the differences in terms of magnitude units (m.u.) are also similar between the other SOF. This suggests that a large part of the SOF term is influenced by radiation pattern effects, and interpretations of these differences in terms of geodynamics or hazard potential should be done very cautiously.

Figure 11

$M_{\mathrm{e}}$ versus $M_{\mathrm{e}}\left(\mathrm{BB}\right)$ and $M_{\mathrm{e}}\left(\mathrm{HF}\right)$, as categorized with SOF.

[Figure omitted. See PDF]

5 Real-time module for SeisComP

The module, derived from me-compute has been integrated to the SeisComP package () and has been part of the GEOFON routine real-time processing since December 2021. The first event for which $M_{\mathrm{e}}$ calculations are available and disseminated via the usual GEOFON services is https://geofon.gfz-potsdam.de/eqinfo/event.php?id=gfz2021xxzt (last access: March 2024​​​​​​​), which occurred on 7 December 2021 10:28:00.3 UTC ($M_{\mathrm{e}}$ 5.7 and $M_{\mathrm{w}}$ 5.5). The scmert add-on is available at https://github.com/SeisComP/scmert (last access: March 2024​​​​​​​).

Figure 12

Screenshot of the SeisComP Origin Locator View (scolv) interactive tool for the $M_{\mathrm{w}}$ 7.7 Republic of Türkiye earthquake that occurred on 6 February 2023, 01:17 UTC along the East Anatolian Fault. The obtained network magnitude value of $M_{\mathrm{e}}$ is 7.8. Stations used are color-coded according to $M_{\mathrm{e}}$ magnitude residuals (top-left frame). Stations in gray are excluded from the network magnitude, do not match the distance range definition, or were trimmed while computing the average magnitude because they are within the $\pm \mathrm{12.5}$ % range. The top-right scatterplot shows $M_{\mathrm{e}}$ residuals by distance (in red those that contributed to actual $M_{\mathrm{e}}$ network magnitude). The topography shown in the map is generated using the ETOPO1 global relief model .

[Figure omitted. See PDF]

The add-on has been configured in GEOFON to trigger the calculation for each origin created by the automatic processing with magnitude $\ge \mathrm{5.5}$ and to compute station magnitudes $M_{\mathrm{e}ij}$ for all stations according to the definition of $M_{\mathrm{e}}$ in the distance 20–98°. The scmert procedure is applied with the settings used by the GEOFON Earthquake Information Service, using stations available in real time from the GEOFON Extended Virtual Network (https://geofon.gfz-potsdam.de/eqinfo/gevn/, last access: March 2024), including station selection and the distribution trimming of 25 %. The workflow for $M_{\mathrm{e}}$ computations is as follows: as soon as an automatically detected event reaches the magnitude threshold, scmert is triggered and starts to compute $M_{\mathrm{e}ij}$ upon receiving data from stations beyond 20°. The process continues until the selected window length (determined by the actual preliminary magnitude) of the last station at 98° is acquired. The first estimate of the magnitude $M_{\mathrm{e}}$ is released shortly after collecting 20 $M_{\mathrm{e}ij}$ estimates from individual station, usually within a few minutes of the earthquake's origin time. SeisComP modules continue to refine the estimate until no further updates are required (this includes manual release at later stages). The computed station magnitudes $M_{\mathrm{e}ij}$ are fully integrated also into the SeisComP Origin Locator View Graphical User Interface (scolv GUI; Fig. ), with station magnitudes and residuals displayed in a dedicated energy magnitude tab.

The energy magnitude values from both modules are compared in Fig. . We used scmert with the same settings as the GEOFON earthquake monitoring service, including station selection and trimming of the distributions. The values are in good agreement, and the best-fit model is $M_{\mathrm{e}}=\mathrm{0.057}+\mathrm{0.987}{M}_{\mathrm{e}}\left(\mathrm{GEO}\right)$, with a standard deviation of 0.118. The average difference computed for magnitudes between 6 and 8 is $-\mathrm{0.028}$.

All values for $M_{\mathrm{e}}$ that have been calculated since the start of the routine processing with scmert can be accessed via the fdsnws-event web service running at GEOFON by specifying $M_{\mathrm{e}}$ as the magnitude type (i.e., https://geofon.gfz-potsdam.de/fdsnws/event/1/query?starttime=2021-12-07&magnitudetype=Me&includeallmagnitudes=true&nodata=404, last access: March 2024). These values are also disseminated to other agencies (e.g., ISC and EMSC) via the usual downstream channels, including a real-time push service.

Figure 13

Comparison between $M_{\mathrm{e}}$ computed in real time by GEOFON, with the scmert add-on for SeisComP ($x$ axis) and offline estimation using me-compute ($y$ axis), considering 153 common events.

[Figure omitted. See PDF]

6 Code and data availability

Code used for computing the energy magnitude is available from

• offline computations in me-compute https://doi.org/10.5880/GFZ.2.6.2023.008 ; and

• real-time computations in SeisComP with scmert https://doi.org/10.5880/GFZ.2.4.2020.003 .

Analyses have been performed in R (https://www.R-project.org/, ), and we used the Generic Mapping Tools (https://www.generic-mapping-tools.org/, last access: March 2024; ) to produce Figs. 2, 5, 6, and 7. The archive, including the energy magnitude catalog (D3 and D6 in Table 1), the complete list of references for the seismic networks analyzed in this article with me-compute, and an example of the configuration files, is available from (10.5880/GFZ.2.6.2023.010).

7 Conclusive remarks

We computed the energy magnitude $M_{\mathrm{e}}$ for 6349 events in the moment magnitude catalog disseminated by GEOFON. When combined with $M_{\mathrm{w}}$, $M_{\mathrm{e}}$ allows a better characterization of the tsunami and shaking potential of an earthquake. The procedure used to compile the dataset, which includes 1 031 396 $M_{\mathrm{e}}$ values for each recording station, is described in detail. Residuals are evaluated using a mixed-effects regression, which partitions the overall residuals into event-specific and station-specific contributions. These random effects are included in the distributed catalog, enabling the computation of $M_{\mathrm{e}}$ for future events and using inter-station residuals as station corrections to reduce the uncertainty in $M_{\mathrm{e}}$. They also enable the assessment of energy magnitude adjustments for specific regions or faulting mechanisms, using interevent residuals and locating propagation anomalies with respect to the global model used to compute Green's functions using the leftover residuals. The methodology employed for computing $M_{\mathrm{e}}$ is suitable for the rapid assessment of $M_{\mathrm{e}}$ . Therefore, it has been implemented as a module for SeisComP, allowing the automatic computation of $M_{\mathrm{e}}$ in real time and keeping the $M_{\mathrm{e}}$ catalog up to date.

Author contributions

DB, AS, and DDG conceptualized the study. RZ developed the Python code used to compile the disseminated catalog. AH developed the add-on for SeisComP. DB developed the quality checks. PE, AH, and AS organized the publication of $M_{\mathrm{e}}$ by GEOFON via web services. All authors participated to the finalization of the article.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We thank all network operators providing data via EIDA–ORFEUS and IRIS, as well as all real-time data providers contributing to the GEOFON virtual network. We thank Kirsten Elger for her help in preparing the dissemination of the products of this study. Finally, we appreciated the valuable feedback and suggestions provided by two anonymous reviewers and by the editor Andrea Rovida.

Financial support

he article processing charges for this open-access publication were covered by the Helmholtz Centre Potsdam – GFZ German Research Centre for Geosciences through the finacial support of the Horizon Europe project Geo-INQUIRE, funded by the European Commission (HORIZON-INFRA-2021-SERV-01, project no. 101058518).The article processing charges for this open-access publication were covered by the Helmholtz Centre Potsdam – GFZ German Research Centre for Geosciences.

Review statement

This paper was edited by Andrea Rovida and reviewed by two anonymous referees.