# Qseek: A data-driven Framework for Automated Earthquake Detection, Localization and Characterization

Marius Isken, Peter Niemz, Jannes Münchmeyer, Pinar Büyükakpinar, Sebastian Heimann, Simone Cesca, Hannes Vasyura-Bathke, Torsten Dahm

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marius Isken, Peter Niemz, Jannes Münchmeyer, Pinar Büyükakpinar, Sebastian Heimann, et al. "Qseek: A data-driven Framework for Automated Earthquake Detection, Localization and Characterization. Seismica, 2025, 4 (1), 10.26443/seismica.v4i1.1283. hal-04976482

HAL Id: hal-04976482

https://cnrs.hal.science/hal-04976482v1

Submitted on 22 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Qseek: A Data-Driven Framework for Automated Earthquake Detection, Localization and Characterization

Marius Paul Isken ID * 1,2, Peter Niemz ID 1,3, Jannes Münchmeyer ID 4, Pınar Büyükakpınar, ID 1, Sebastian Heimann ID 2, Simone Cesca ID 1, Hannes Vasyura-Bathke ID 1, Torsten Dahm ID 1,2

¹Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences, Potsdam, Germany, ²University of Potsdam, Germany, Institute of Geosciences, Potsdam, Germany, ³University of Utah, Seismograph Stations, Salt Lake City, UT, USA, ⁴Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTerre, Grenoble, France

**Author contributions:** Conceptualization: Marius P. Isken, Sebastian Heimann. Methodology: Marius P. Isken, Sebastian Heimann, Simone Cesca. Software: Marius P. Isken, Sebastian Heimann. Validation: Hannes Vasyura-Bathke, Simone Cesca, Peter Niemz, Pınar Büyükakpınar. Formal Analysis: Jannes Münchmeyer, Sebastian Heimann. Writing - original draft: Maius Paul Isken. Writing - Review & Editing: Marius Paul Isken, Torsten Dahm, Simone Cesca, Peter Niemz, Jannes Münchmeyer, Pınar Büyükakpınar. Visualization: Marius P. Isken, Peter Niemz. Supervision: Torsten Dahm, Simone Cesca. Funding acquisition: Torsten Dahm.

**Abstract** We introduce a data-driven method and software for detecting and locating earthquakes in large seismic datasets. By combining seismic phase arrival annotations, delivered by neural network phase pickers, and waveform stacking with an adaptive octree search, we can automatically detect and locate seismic events even in noise-dominant seismic data. The resolution of the search volume is iteratively refined toward the seismic source location; this strategy facilitates an efficient, fast, and accurate search. We present a user-friendly and high-performance open-source software framework based on established frameworks, featuring event detection in layered 1D and complex 3D velocity models and event feature extraction capabilities, such as moment and local magnitude calculation from peak ground motions. We incorporated station-specific corrections and source-specific station terms into the search to enhance the location accuracy. We demonstrate and validate our approach by extracting extensive earthquake catalogs from large seismic datasets in different regions and geological settings: (1) Reykjavanes Peninsula, Iceland; (2) Eifel volcanic region, Germany; and (3) Utah FORGE, USA. We capture seismic events from tectonic activity, volcanic swarms, and induced microseismic activity with magnitudes ranging from -1 to 5. Such precise and complete earthquake catalogs contribute to the interpretation and understanding of otherwise hidden subsurface processes.

**Non-technical summary** We present a new method and open-source software for automatically detecting and localizing earthquakes and microseismicity in large seismic datasets recorded by seismometers. The technique can detect natural earthquakes, such as tectonic faulting events, volcano-tectonic swarm activity, and induced seismicity from well operation in oil and gas or geothermal exploitation. Our method uses machine learning and an adaptive focusing mechanism to efficiently search through large amounts of continuous seismic data and detect and locate earthquake events. We present an open-source software written in Python, qseek, for this purpose. To test and showcase the method, we look at seismic data from (1) Iceland, where a volcano-tectonic sequence was recorded, (2) a large dataset from the Eifel Volcanic Region, Germany, and (3) microseismicity recorded at Utah FORGE, USA, a field-scale geothermal lab. By detecting and locating the small seismicity underground, we make otherwise hidden processes within the Earth visible. Our robust and easy-to-use method contributes to understanding natural seismic activity and manmade seismicity.

## 1 Introduction

The detection and localization of seismic events provide information about deformation processes and rock volumes under high elastic stress in the Earth's brittle crust. The (micro)seismicity often stems from tectonic and volcanic processes but can also occur in response to human activities such as fluid extraction or injection (Ellsworth, 2013; Vasyura-Bathke et al., 2023; Niemz et al., 2024). By accurately mapping seismic activity in space and time, we gain deeper insights into

these dynamic processes, subsurface structures, and the state of stress in the Earth (e.g. Cesca et al., 2020; del Fresno et al., 2023; Greenfield et al., 2022; Wilding et al., 2023).

The advent of dense permanent and large-N temporary seismic monitoring networks has markedly enhanced our observational capabilities, resulting in a growing influx of seismic data, capturing a spectrum of seismic events (e.g. Dougherty et al., 2019; Obermann et al., 2022; Dahm and The Eifel Large-N team, 2023). Among these are faint microseismic events, often obscured by seismic noise and characterized by a low

Production Editor:

Carmine Galasso

Handling Editor:

Carmine Galasso

Copy & Layout Editor:

Hannah F. Mark

Received:

February 29, 2024

Accepted:

January 21, 2025

Published:

February 19, 2025

<footnote>*Corresponding author: marius.isken@gfz.de</footnote>
signal-to-noise ratio (SNR). Traditional seismic detection and localization methods, as well as recently presented phase associators for large numbers of picks (e.g. Zhu et al., 2022; Münchmeyer, 2024) often struggle with noise-dominant micro-seismic events or massive concurrent occurrences (e.g. Grigoli et al., 2013; Li et al., 2020).

To address the challenge of detecting and locating weak seismic events, waveform-stacking methods have emerged that enable robust event detection and localization through energy stacking (Gharti et al., 2010; Cesca and Grigoli, 2015; Grigoli et al., 2013). Several adapted versions of this backprojection method have been suggested and applied to various purposes, such as observing natural tectonic sequences and monitoring of volcano-tectonic (VT) swarms or induced seismic activity (Greenfield et al., 2020; Cesca et al., 2020; del Fresno et al., 2023; Wilding et al., 2023; Flóvenz et al., 2022). The key benefit of these data-driven methods is their ability to produce reliable localizations (e.g. Kao and Shan, 2004; Drew et al., 2013; Poiata et al., 2016; Grigoli et al., 2013; Winder et al., 2021). Notable established software frameworks for stacking and migration are QuakeMigrate (Winder et al., 2021), LOKI (Grigoli et al., 2013), its companion MALMI (Shi et al., 2022) and BPMF (Beaucé et al., 2023). Some older frameworks suffer performance issues when confronted with large datasets and extensive search volumes. These scaling issues are one of the main motivators for developing the method presented here.

We introduce a fast, fully automatic method and a scalable, high-performance software framework driven by waveform data and machine learning (ML) to detect, locate, and characterize seismic events in large, continuous seismic datasets. The method is based on stacking and migration, which capitalizes on the coherency of seismic arrivals to reveal obscured low-magnitude events. Seismic energy arrivals at individual stations are annotated by neural networks (NN) trained to identify phase arrivals in seismograms (Zhu and Beroza, 2018; Münchmeyer et al., 2022). A volume grid search is used to stack and correlate individual arrivals of seismic energy based on a travel time model. To further enhance the efficiency and convergence of the search, we replace the regular grid with a hierarchical octree. The octree can be adaptively refined towards the seismic source region, resulting in an efficient search and more precise event locations while reducing computational effort. Subsequently, station-specific correction terms (SST) can be obtained by extracting time delays between calculated and observed phase arrivals. Source-specific station terms (SSST) can also be computed to account for subsurface heterogeneities. These SST and 3D SSST station corrections can subsequently be applied in a second detection and localization search to augment the travel time model, resulting in more precise localizations and more event detections (Nooshiri et al., 2016; Grigoli et al., 2016; Lomax and Savvaidis, 2022).

Along with the method, we present the modular open-source, high-performance software framework qseek, written in the Python and C programming languages for handling and analyzing large continuous

seismic waveform datasets. It features forward-modeling of seismic travel times in layered 1D and complex 3D Earth velocity models. It utilizes SeisBench (Woollam et al., 2022), which provides a uniform interface to different pre-trained neural-network-based phase annotators. Moreover, the software framework enables the automatic extraction of event features, such as local and moment magnitudes from peak amplitudes, and station corrections in the form of SST and SSST. The framework is flexible to different use cases and seismological applications.

To demonstrate the effectiveness of this approach, we present the earthquake catalogs extracted from three distinct seismic datasets: (1) a volcano-tectonic dataset from the Reykjanes Peninsula, Iceland, featuring VT swarm seismicity during a period of episodic magmatic unrest in 2020, (2) a large-N seismic dataset from the Eifel Volcanic Region, Germany, where both tectonic earthquakes and seismic swarm activity–evidence of the volcanic activity in the region–are recorded, and (3) an example from Utah FORGE, covering a case of induced microseismic activity caused by the circulation of fluids within the subsurface, monitored by a small seismic network.

This work focuses on method development and leaves the detailed geological interpretation of the obtained seismic catalogs to future studies. The reader is directed to further publications that provide an in-depth analysis of the seismicity catalogs presented here.

## 2 Method

The earthquake detection and localization method and software framework consist of six key components: (1) the waveform image function annotating phase arrivals, (2) the seismic travel time model, (3) stacking and migration, (4) the octree volume, (5) station corrections, and (6) event feature extraction. Some components will receive more significant focus, while others are only briefly introduced with references to existing foundational works.

### 2.1 Waveform Image Function: Machine Learning First Arrival Annotation

The recorded waveforms are transformed into non-negative image functions (or characteristic functions) to stack the detected P and S wave arrivals constructively. These image functions (IF) require sensitivity to P and S wave onset energy to image the first arrivals (i.e., phase picks). This can be achieved by analytical characteristic functions (CF), typically variants of the STA/LTA function or the waveform envelope (Cesca and Grigoli, 2015; Grigoli et al., 2014; Lomax et al., 2012). In this study, we replace the conventional analytical image functions with a neural network (NN) trained to identify P and S wave onsets. These established machine learning annotators have higher specificity and sensitivity than simple analytical CFs (Shi et al., 2022).

The emergence of user-friendly machine learning frameworks such as TensorFlow and PyTorch has facilitated the design and efficient training of neural networks.
works for seismic phase detection and arrival picking (Zhu and Beroza, 2018; Paszke et al., 2019). These networks are trained on large volumes of labeled datasets from global seismological entities and observatories. They can detect phase arrivals in single station seismograms D with high sensitivity and pick precision (Zhu and Beroza, 2018; Mousavi et al., 2020; Münchmeyer et al., 2022). We utilize the neural networks annotated output (NN) as an image function I for P and S wave arrivals (Fig. 1, Shi et al., 2022). Seis-Bench (Woollam et al., 2022) is used as a uniform and performant software framework to access various pre-trained neural networks for seismic phase annotation, such as PhaseNet, EQTransformer, and GPD (Zhu and Beroza, 2018; Mousavi et al., 2020; Ross et al., 2018). Münchmeyer et al. (2022) have comprehensively evaluated and compared various NN architectures and training datasets.

The image function I(t) implemented in qseek is given by

$$I^{\theta}(t) = NN^{\theta}(D(t)) \quad \theta \in \{P, S\} \qquad (1)$$

where NN is the annotating neural network, θ is the seismic phase and D(t) is the recorded seismogram.

<figure><img src="image_3.png" /><figcaption>**Figure 1** (a) Single station seismic waveforms from a local earthquake recorded on Reykjanes Peninsula, Iceland. The three-component seismogram D(t) is input to the neural network image function. (b) The neural network's (Seis-Bench/PhaseNet) annotations of P and S phase arrivals, i.e., the image function I(t) and peak-detections (pick times; dashed lines). The annotated output has a Gaussian shape with a half-width of 0.2 seconds; the maximum value represents the confidence of the annotated pick. The analytical P- and S-STA/LTA characteristic functions for the vertical (window ratio 1:5) and horizontal traces (10:30) are shown in the background. Following the calculation of Grigoli et al. (2014).</figcaption></figure>

The image function produced by the neural network

presents the advantage that phase arrival annotations are represented as peaks with a Gaussian-like shape with a standard deviation of 0.2 seconds (e.g. Münchmeyer et al., 2022). The peak value provides a measure of phase identification certainty, while the shape provides certainty in time. The peak value max(I(t)) is a qualitative indicator and can be interpreted as a confidence measure p, ranging from 0 to 1 (Fig. 1). The NN image function's Gaussian-like shape and maximum value are less ambiguous than the erratic STA/LTA trigger (Allen, 1982). Empirical observations suggest that NN residuals are symmetric around the actual pick, while STA/LTA characteristic functions usually have higher values after the true arrival. Notably, neural networks are commonly trained with 100 Hz input data, resulting in 100 Hz annotation functions. (Zhu and Beroza, 2018; Mousavi et al., 2020). To improve the annotation sensitivity for, e.g., microseismic events, the self-similarity character of earthquake signals can be exploited by stretching and compressing the waveforms, thereby shifting (or rescaling) the effective frequency content of the input waveform data D(t) with a constant factor Sᵢ (Shi et al., 2022).

## 2.2 Seismic Travel Time Model

Calculation of seismic travel times τ<sub>mod</sub><sup>θ</sup> from potential event locations x to the receiving station s is required to accurately shift phase arrival times for stacking and subsequently detect and locate event hypocenters (e.g. Drew et al., 2013). The calculation of travel times relies on an apriori velocity model. This model can be approximated as a layered 1D or a 3D model and inferred from wave travel time analysis, seismic tomography, and stratigraphic information. In the following, we limit the explanations to direct P and S phases since the image functions (eq. 1) identify P and S phase arrival times. The wave travel time τ from an event with origin time T₀ to a station s can be described as:

$$\tau_s^\theta(\mathbf{x}) = A_s^\theta(\mathbf{x}) - T_0(\mathbf{x}) \quad \theta \in \{P, S\}, \qquad (2)$$

where A_s is the absolute arrival time of the seismic phase θ (P and S direct phases) at the station s (Stein and Wysession, 2002).

### 2.2.1 1D Layered Velocity Model

While 1D layered velocity models simplify the complex 3D structure of the Earth, they provide valuable insights into wave propagation and wave physics. This simplification is a first-order approximation of the subsurface, which is computationally effective and practical, mainly when high-resolution 3D velocity models are unavailable. For calculating phase travel times τ in 1D layered models, we utilize Pyroko's cake ray tracer (Heimann et al., 2017) based on classical ray theory. Ray tracers are essential for comprehending how seismic waves spread, refract, and reflect within the Earth. They are also able to forward model complex ray geometries and phases.
### 2.2.2 3D Velocity Model: Fast-Marching Eikonal Solver

In seismology, travel times and ray paths for direct waves in 3D heterogeneous media can be calculated with the Eikonal equation (Hamilton, 1828; Dimri et al., 2003) to achieve more precise travel time calculations. We use the efficient numerical fast marching method to obtain wave travel times for a search region on an evenly spaced grid (Sethian, 1996). Our method implementation can process and read NonLinLoc velocity grids (Lomax et al., 2012). We calculate inverted travel times from the receiver towards all grid nodes to efficiently compute and store travel times between all seismic stations and all grid nodes (Fig. 2). This allows the creation of interpolated travel time volumes and the effective caching of node travel times without sacrificing computational performance.

<figure><img src="image_4.png" /><figcaption>**Figure 2** Conceptual 2D visualization for seismic travel time calculation in heterogeneous velocity media using the fast-marching method for the Eikonal solution. Inverted travel times from the receiving station at the surface (indicated by a yellow triangle) towards all grid nodes in the subsurface are calculated, simultaneously resulting in station-specific travel times for all potential source locations. The contour lines indicate the travel time.</figcaption></figure>

## 2.3 Search by Stacking and Migration

Stacking and migration is an established beamforming method to detect and locate seismic sources within a defined volume. By time-shifting positive image functions according to a calculated travel time $\tau_s^\theta$ between the grid node location x in the subsurface and the observation location s, we effectively backproject the seismic energy towards the potential source location x. In that way, we time-shift and stack image functions from many observations (stations) onto each grid node within the search volume, yielding the quantity called the semblance $E$. Then the point of maximum semblance $E(x, t)$ (also known as coalescence, coherence or beam; Drew et al., 2013; Grigoli et al., 2013; Beaucé et al., 2023) within the search volume indicates the most likely source location at time $t$. This method involves a grid search for the maximum semblance across all grid points inside the source volume. When energy enters the search volume

from the outside, maximum semblance will be located at or inside the volume's border. To counteract this, we incorporate an absorbing boundary. Events within this boundary are ignored or can be filtered out for analysis. Inaccurate seismic velocity models result in a less constructive stack, leading to fewer detections and the deterioration of localization accuracy. A comprehensive description of the waveform stacking method can be found in Drew et al. (2013) or Grigoli et al. (2014). The stacking and migration concept for a single node x in a 2D grid is shown in Figure 3. The extension to a 3D grid follows the same concept. The semblance $E$ at a grid location x, time $t$ for a seismic phase $\theta \in \{P, S\}$ is given by Grigoli et al. (2016)

$$E_{\theta}(\mathbf{x}, t) = \sum_{s}^{S} I_{s}^{\theta}(t + \tau_{s}^{\theta}(\mathbf{x})) \quad (3)$$

where $\tau_s^\theta$ is the theoretical travel time of a seismic phase at seismic station s in a set of stations S.

<figure><img src="image_5.png" /><figcaption>**Figure 3** Concept of stacking and migration grid search method outlined in 2D: The positive seismic trace image functions $I$ are shifted by calculated travel-times ($\tau_s^\theta$) and stacked for each node x in the grid. For each grid point, we obtain a stacked time series of semblances. This semblance stack reaches its maximum at the most probable hypocenter location at event time $t$ of the seismic event, sketched by the star.</figcaption></figure>

The semblance functions $E_\theta$ for P and S phase arrivals are combined to a semblance function $\hat{E}$ using the power mean (eq. 4), where $\eta$ is a tune-able parameter that shapes the stacked semblance: Values $\eta > 1$ will sharpen the individual and combined semblance and result in higher location accuracy; and $\eta < 1$ will smooth the stacked image function, which can be of advantage when the velocity model is insufficient to stack the different image's phases. Typical values for $\eta$ range from 1 to 1.5.
The power mean is defined as

$$ \hat{E}(E_{\theta}, \eta) = \left( \frac{\sum_{\theta} w_{\theta} E_{\theta}^{\eta}}{\sum_{\theta} w_{\theta}} \right)^{1/\eta} \quad (4) $$

where $w_\theta$ is the user-defined weight for phase $\theta$. Commonly, P and S phase weight is 1 for well-posed datasets. Other works use the general coherence function to combine the P and S phase semblance functions (e.g. Grigoli et al., 2016), but cannot control the shape of the input image functions.

It is inherent to the location search that stations near the hypocenter provide the most precise information about the event's location. Furthermore, as the distance between stations and the event grows, the ray path complexity rises, leading to more significant travel time modeling errors with distance. We propose a station-distance weighting, $w_s$, to offset these effects based on the node-station distance:

$$ w_s(\mathbf{x}) = \frac{1 - W}{1 + \left(\frac{r(s, \mathbf{x})}{R}\right)^3} + W \quad (5) $$

where $\mathbf{x}$ is the node location, $s$ the station, $r$ the 3D distance from the node to the station, $R$ is the radius of the spatial decay function, usually 10 - 50 km (depending on the search volume and station density) and $W$ is a water level, ensuring that stations outside the cubic decay function are included in the search.

## 2.4 Octree Search

The semblance $\hat{E}(\mathbf{x}, t)$ is defined for every potential grid point $\mathbf{x}$ and time $t$. Local maxima detect events in the semblance stack. To identify these maxima, we follow a two-step procedure. First, we reduce the semblance function to the peak semblance $\bar{E}(t)$ at each time step:

$$ \bar{E}(t) = \max_{\mathbf{x}} \left( \hat{E}(\mathbf{x}, t) \right) \quad (6) $$

We keep track of the location $\mathbf{x}$, maximizing the semblance at each time step $t$. We then identify peaks in the semblance function $\bar{E}(t)$ using a simple trigger criterion. This reduction might lead to missed detections if two events occur with very short interevent time (< 1 s), for example, during borehole injection or very dense seismic swarm activity.

To perform the spatial reduction of the semblance function, we employ a hierarchical octree data structure (Meagher, 1982). The octree has coarse resolution in regions of low semblance and can iteratively adapt to high resolution in areas of high semblance. For a single time step $t$, we initialize our root octree $O$ with a coarse, regular grid covering the entire search volume (Fig. 4). We call each cubical grid cell a node $n$ with edge length $l_n$ and define the semblance $\hat{E}(n, t)$ as the stacked semblance trace at the center point of the node $n$ at time $t$. We then apply an iterative procedure for refining the octree to find the node with the maximum semblance:

**Step 1** Calculate semblance $\hat{E}(n, t)$ for all nodes in the octree (caching can be applied)

**Step 2** If for no node $n \in O$ the predefined semblance threshold $E_0$ is exceeded, i.e. $\forall n \in O: \hat{E}(n, t) \le E_0$, terminate the octree search.

**Step 3** Select nodes: (1) with the maximum semblance $n_m$ and (2) the node with the maximum semblance density $n_v$ per volume, and (3) their neighbor nodes $n_n$ with contacting face, edge, or vertex. The semblance density per volume $\hat{E}^v$ is defined as:

$$ \hat{E}_n^v = \frac{\hat{E}_n}{l_n^{3}} \quad (7) $$

Each node is split into 8 equally sized child nodes that are added to the tree (Fig. 4). A leaf node cannot be split and is removed from the selection. This strategy guarantees an efficient exploration of the semblance space, and the nodes $n_m$ and $n_v$ will eventually converge at the hypocenter.

**Step 4** If the maximum semblance $\bar{E}(t)$ falls into a leaf node $n_m$ and the selected nodes cannot be split, the hypocenter has been found, and the search can be terminated. Otherwise, go to step 1 and continue searching for refined octree.

**Step 5 (Optional)** Perform a neighborhood maximum search around the maximum semblance node $n_m$ by interpolating $\bar{E}$ using radial basis functions. This will refine the hypocenter location for larger leaf node sizes.

The maximum semblance $\bar{E}(t)$ can now be determined in the same way as before. Note that this optimization might lead to a higher total number of explored nodes but not a change in the maximum semblance function. It is computationally efficient because larger batches are examined, and it enables vectorization. The octree approach facilitates efficient global exploration of expansive volumes and datasets and adaptively refines its resolution towards the seismic event hypocenter by iteratively focusing the search towards regions with pronounced semblance signatures (Fig. 5). Leveraging this strategy accelerates the resource-intensive migration and stacking method, thereby optimizing the scanning process for seismic events within large datasets and large volumes while enhancing the resolution near the hypocenter, leading to more accurate event locations. The effective octree strategy is also employed by NonLinLoc for earthquake localization (Lomax et al., 2000).

For regional analysis, the root level nodes have a typical edge length of 1-2 km with 4-6 refinement levels $l$ exposing bottom nodes with an edge length of 62.5 m. The finest resolution of nodes is limited by the semblance sampling frequency and the minimum resolvable wavelength at the source depth (depending on the seismic velocity at depth). Generally, a higher semblance sampling frequency allows for the resolution of smaller octree nodes and leads to a more accurate semblance stack and better hypocenter locations. On the other hand, root nodes ($l = 0$) that are too large will suffer from aliasing effects in the semblance stack, and peaks in the IF will not stack constructively. Effective
anti-aliasing measures before stacking are reducing the image function with a maximum filter or the convolution with a boxcar function.

The efficacy of a search is dependent upon the semblance threshold $E_0$ to detect events on the octree. A too-small threshold will result in many false positives, while a too-large threshold will result in missed events. The number of total stations, station density, parametrization of the station-distance weighting $w_s$, and ML phase annotation model influence the $\hat{E}(n, t)$ value. It is generally advisable to begin with a lower $E_0$ value close to zero and filter out erroneous detections.

<figure><img src="image_6.png" /><figcaption>**Figure 4** Refinement of an octree: one cubic root node (octant) is split into eight smaller, equally-sized nodes. The illustration displays a two-level refinement of a single node from left to right.</figcaption></figure>

## 2.5 Location Uncertainty

Uncertainties of the detected event locations can be determined in different ways: (1) Perturbing the image function and velocity model and relocating the event many times (Grigoli et al., 2013). This method is similar to Bayesian bootstrapping and computationally costly. (2) The temporal-spatial semblance function can be fitted with a multivariate Gaussian distribution (Drew et al., 2013). (3) Transforming the semblance into a probability density function (Poiata et al., 2016).

We propose extracting the event's location uncertainty from the spatial semblance distribution of surrounding octree nodes: Nodes are selected within a pre-defined semblance percentile, usually set at 2%. The spatial extent of these selected nodes defines the location error in all three dimensions, thus providing an orthonormal triaxial error ellipsoid for the location. See Figure 5 for a 2D projection of the node semblance distribution around the located hypocenter. This method is simple yet a robust qualitative measure for the uncertainty of the location but has shortcomings in non-ellipsoidal, non-orthonormal, and multimodal semblance distributions. An example of the uncertainty procedure is demonstrated for seismicity on the Reykjanes Peninsula, Iceland, and Eifel Volcanic Region, Germany (see sections 3.1 and ??).

## 2.6 Station Corrections

Simple 1D and 3D seismic velocity models are often不适 to represent geological reality. This leads to inaccuracies between the forward-modeled theoretical $\tau_s^{\theta}$ and observed arrival times $t_{\text{obs}}^{\theta}$ of P and S first arrivals (Billings et al., 1994). These inaccuracies reduce the accurate and constructive stacking of the image functions ($E$), leading to inaccurate event hypocenter locations (Grigoli et al., 2016; Lomax and Savvaidis, 2022).

To address these discrepancies, station travel time corrections $\delta\tau_s^{\theta}$ for first arrivals (P and S) can be statistically derived from a catalog of observed travel time delays. In the following, we will describe two variants of the station-specific delay times (or station corrections): simple scalar station-specific terms (SST) and 3D source-specific station terms (SSST).

In a subsequent scan of the seismic dataset, the extracted station corrections can be added to the calculated travel times $\tau_s^{\theta}$ to enhance the precision of the seismic event locations. This augmentation increases the location accuracy and the semblance value, subsequently increasing the number of event detections without increasing the number of false detections, as demonstrated by Grigoli et al. (2016).

### 2.6.1 Station Specific Terms (SST)

The modeled travel time $\tau_s^{\theta}$ and the phase arrival time $t_s^{\theta}$ obtained from the neural network annotation are known. Thus, the average phase-specific residual delay time $\delta\tau_s^{\theta}$ at station $s$ can be statistically extracted from a catalog of located events (Fig. 6):

$$ \delta\tau_s^{\theta} = \tau_s^{\theta} - t_s^{\theta} \qquad (8) $$

The semblance value is proportional to the location confidence or the accuracy of the event. The station-specific correction time is then calculated by the weighted mean or weighted median of all observed travel time delays $\delta\tau_s^{\theta}$:

$$ \delta(\tau_s^{\theta}) = \text{stat}(\delta\tau_{s,1}^{\theta}, \dots, \delta\tau_{s,n}^{\theta}, w_i^p), \quad \text{where } \text{stat} \in \{\text{mean, median}\} \qquad (9) $$

and

$$ w_i^p = p_i \bar{E}_i \qquad (10) $$

where the pick weight $w_i^p$ is defined by the confidence of the arrival pick $p_i$ and semblance $E_i$ of the event detection $i$.

Simple scalar SST station corrections are adequate for compensating common ray path delays beneath the station's location. However, multimodal station delay statistics hint towards more complex 3D heterogeneities of the subsurface (Lin and Shearer, 2005).

### 2.6.2 Source-Specific Station Terms (SSST)

We can calculate source-specific station terms to account for the natural 3D variability of seismic velocities in the Earth. These are station-specific 3D volumes of correction times (Nooshiri et al., 2016; Lomax and Savvaidis, 2022). For the volume, we re-use the octree grid fixed at a defined level $l$ (e.g., $l = 0$); this effectively represents a regular grid. The delay time for a location can be calculated as described in 2.6.1 with an additional node-event distance weighting for each observed event. This distance weighting $w_e(r)$ can be as described as:

$$ w_e(r) = \frac{1}{1 + \left(\frac{r}{R_n}\right)^b} \qquad (11) $$
<figure><img src="image_7.png" /><figcaption>**Figure 5** Horizontal and vertical surface reductions (projections) of the 3D octree showing the first three refinement levels during detection and adaptively focusing onto the seismic source region (left to right). The light color indicates high semblance, and the dark indicates low semblance. At level 0 (left) the edge length of the nodes is 1 km, these are split into smaller nodes of 500 m and 250 m in length at refinement levels 1 and 2. The 3D node count increases from 5977 to 6812 nodes for the levels shown. The shown earthquake detection is for a regional event detected within the Eifel large-N network. Not illustrated here are higher refinement levels 3 (125 m) and 4 (62.5 m).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6** Distribution of SST delay times (red negative, blue positive) at station ELB annotated with different statistical measures of the distribution: The weighted (wt.) mean and median are calculated by using weights based on the pick confidence multiplied with the event's semblance.</figcaption></figure>

where r is the distance between node centroid and event hypocenter, and $R_n$ is the size of the sphere enveloping the node. The spatial decay exponent b is typically 3. The distance weighting can be combined with the scaled pick weight $w_i^p$ from eq. 10. Then we find the radius $R_n$ around each node to encompass picks with integrated pick confidence $P_n$ greater than threshold $P_T$:

$$P_n = \sum_{i}^{K} p_n$$
$$\text{where } K = \min\{k \in \{1, 2, \dots, n\} : S_k > P_T\} \quad (12)$$

This approach is similar to Richards-Dinger and Shearer (2000) but recognizes the phase pick confidences to find a sufficiently large $R_n$ for each node. For our analysis, we chose $P_T = 5$. The obtained delay time volume (Figs. 7, S3) has an identical resolution as the octree at refinement level l. The station and phase-specific delay 3D grid of phase delay times can then be used to interpolate and extrapolate station correction times for arbitrary source locations.

## 2.7 Event Feature Extraction

After detecting and localizing seismic events, we can extract event features from the seismic waveform data. These may include event magnitudes or ground motion observations that can be input for shake maps or site amplification maps.

### 2.7.1 Local Magnitudes

The calculation of earthquake local magnitude $M_L$ from measured peak ground motions follows different regional and local magnitude definitions (Bormann and Dewey, 2012). This complexity of local magnitudes is
<figure><img src="image_9.png" /><figcaption>**Figure 7** SSST delay times of the P first arrivals for a single seismic station SI.VOS on Reykjanes Peninsula, Iceland. Slices are from the full volume (Fig. S3) reduced to (a) a horizontal 2D surface and (b) a vertical cross-section covering the octree volume. The SSST residual times are obtained from weighted travel time residuals calculated on the octree at level $l = 1$ (1 km grid spacing) and $P_T$ of 5. Travel time residuals shown here are calculated from 50k+ events on the Reykjanes Peninsula (black dots).</figcaption></figure>

inherent in the local attenuation model and the model's definition of peak ground motions $A_P$ measurement (i.e., horizontal or vertical components). Most attenuation models require the simulation of a Wood-Anderson seismometer; the accepted numerical representation of this response has historically changed (Hutton and Boore, 1987; Uhrhammer and Collins, 1990; Bormann and Dewey, 2012), while other attenuation models require the measured particle motion (e.g. Horálek et al., 2000). The qseek software framework has incorporated various local magnitude models for estimating event $M_L$ in different regions. The framework estimates the individual station magnitudes $M_L^S$ together with pre-event noise $D_N$, noise standard deviations $\sigma(D_N)$, noise amplitudes $A_N$ and peak amplitudes $A_P$. Valid station peak amplitudes $A_P$ are selected by $A_P > 2\sigma(D_N)$. Clipped waveforms (digital and analog) are filtered out. Magnitude errors can then be calculated by adding and subtracting the noise amplitude $\pm A_N$ to the peak amplitude estimation $A_P$ and recalculating the station's local magnitude. The network's local magnitude $M_L$ is calculated from the median of estimated station magnitudes $M_L^S$ (Bormann and Dewey, 2012). Individual station magnitude corrections, however, are not estimated, and the local magnitudes calculated here are less reliable than those derived from a calibrated seismic network.

### 2.7.2 Moment Magnitudes

A data-driven method for moment magnitude estimation from modeled peak ground motions is proposed by Dahm et al. (2024). It retains traditional magnitude determination from peak amplitudes while addressing

limitations like rigid frequency filters and responses, signal saturation, or source complexity.

The relationship of moment magnitude $M_W$ and observed peak ground motion is given by Dahm et al. (2024):

$$M_W = \text{median} \left( \log_{10} \left( \frac{u_s}{\hat{u}^R} \right) + M_W^R \right). \quad (13)$$

where $u_s$ is the measured peak ground displacement (or velocity) at station s and $\hat{u}^R$ is the synthetic peak ground displacement and $M_W^R$ a selected reference magnitude to calculate synthetic Green's functions.

The statistical attenuation functions $u(r)$ are forward-modeled through synthetic waveforms provided by Green's function (GF) databases (Heimann et al., 2019). The velocity model underlying the GFs was chosen as the same model used for detecting and localizing the event. This effectively harmonizes the detection, localization, and magnitude estimation procedure.

## 2.8 Software Implementation

The open-source qseek software framework is implemented in Python and C programming languages. The exposed API interface is in Python, while heavy-lifting algorithms (i.e., stacking and migration) are implemented in C. This computationally expensive function is implemented in Pyrocko's parstack module and benefits from massive parallelization on the CPU register level (Single instruction, multiple data; SIMD). Seismic waveform handling is expedited by Pyrocko's squirrel, enabling fast and asynchronous access to large seismic datasets. Node travel time calculations are cached in last-recently-used (LRU) caches. New entries are added to the cache and discarded when not accessed in a long time. This strategy alleviates costly and repetitive forward modeling. An overall modular architecture and structured object-oriented principles organize the code base and allow for the integration of plugins across the entire processing pipeline. This facilitates the integration of other image functions, e.g., analytical phase pickers, such as FilterPicker (Lomax et al., 2012; Bagagli, 2022), or the integration of non-cartesian octrees. To maximize performance and hardware utilization (Disk I/O-RAM-CPU-GPU), we use asynchronous programming principles facilitated by asyncio from the Python standard library.

All data models adhere to semantic schemas and strict typing provided by the Pydantic framework. This enables serialization and de- serialization to and from standard JSON format. The formalization of data exchange formats and semantics, in turn, allows easy handling and examining of large earthquake catalogs. Moreover, it establishes the foundation for data exchange across different software modules (e.g., a web interface or databases). The interplay between different detection and localization workflow components is conceptualized in Figure 8.

### 2.8.1 Waveform Pre-processing

Various frequency filters, such as low-pass, high-pass, and band-pass, can enhance the waveform data's SNR
<figure><img src="image_10.png" /><figcaption>**Figure 8** Components of the presented method and the software framework *qseek* for earthquake detection, localization, and characterization. The components of the method are grouped in blocks (i.e., data, travel time model, stacking, and migration). The calculation and application of station corrections and the extraction of event features are optional modules in the processing workflow.</figcaption></figure>

prior to its annotation by the imaging function. Optionally, the SNR can be improved by applying the Deep-Denoiser neural network (Zhu et al., 2019), which can significantly improve data quality and subsequent NN phase annotation (Münchmeyer et al., 2022). However, pre-processing can be computationally demanding and degrade large datasets' overall search runtime performance.

## 3 Applications

We demonstrate the method and software framework on three seismic datasets containing different scales of seismicity and representing different geological settings: (1) Volcano-tectonic seismic data recorded on Reykjanes Peninsula (RP), Iceland, in 2020, where non-effusive volcanic activity produced episodic swarm seismicity. (2) Seismological data from the Eifel large-N experiment in the Eifel Volcanic Region (EVR), Germany. We detect and locate seismic activity that comprises tectonic earthquakes, small-scale seismic swarm activity, and anthropogenic events. (3) Data from Utah FORGE, a laboratory for enhanced geothermal systems, where we detect and locate microseismicity induced during a circulation experiment in 2023.

In this study, we do not provide interpretations of the located seismicity. For a detailed analysis of the seismicity and source mechanisms observed on RP, we refer to *Büyükakpınar et al.* (2024), for the seismicity in the EVR to *Dahm and The Eifel Large-N team* (2023), and for an analysis of the microseismicity at Utah FORGE to *Niemz et al.* (2024).

### 3.1 Reykjanes Peninsula, Iceland

The Reykjanes Peninsula (RP) in southwestern Iceland represents the intersection of the Reykjanes mid-oceanic Ridge with the island and forms part of an active oblique rift system. Since late 2019, intensified tectonic-magmatic activity has been observed, including earthquake swarms at varying crustal depths, ground deformation, and volcanic fissure eruptions from 2021 to 2024 (Sigmundsson et al., 2024; Flóvenz et al., 2022). These phenomena are linked to fluid migration and magmatic intrusion, activating dike and sill complexes within the crust and upper mantle (Flóvenz et al., 2022;

Ducrocq et al., 2024; *Büyükakpınar et al.*, 2024).

**Analysis** We analyzed seismic datasets comprising 27 stations across the RP and its eastern surroundings for benchmarking and validation. The data span from January 1 to September 1, 2020 (244 days). The earthquake catalog with applied SSST corrections is presented in Figure 9. The catalogs without station corrections and simple SST corrections are shown in supplement Figures S6 and S7. Seismicity maps, including uncertainties, are shown in Figures S8-S12, along with a statistical analysis of uncertainties and station correction terms (Fig. S11). The average horizontal and vertical uncertainties are 233 m and 480 m, respectively (Tab. S2). We compare our results to a curated catalog by the Iceland Meteorological Office (IMO), which is based on a comparable seismic network (Fig. S5). A 3D animation of the seismicity is linked in the supplement.

### 3.2 Eifel Large-N Experiment, Germany

The Eifel Volcanic Region in the Rhenish Massif, Germany, comprises over 350 volcanic centers, including basanitic scoria cones, maars, and phonolitic centers (Mertes, 1983). The last eruptions include the Laacher See Maar (12,800 years ago) (Reinig et al., 2021) and the Ulmener Maar (10,800 years ago) (Zolitschka et al., 1995). The region features active fault systems, such as the SE-NW Ochtendunger fault (ODFZ) with $M < 4$ seismicity and the fossil SW-NE Variscan faults (Hinzen, 2003). Evidence of volcanic activity includes mantle CO$_2$ degassing (Gal et al., 2011; Defourny et al., 2022), upper-crustal swarm earthquakes near older volcanic centers (Hensch et al., 2019), and deep low-frequency swarms at 45 km depth (Dahm et al., 2020b).

In 2022, the GFZ Potsdam and collaborating institutions deployed a temporary large-N seismic network in the German Eifel region focused around the Laacher See volcano (in preparation). The network spans 180 km EW by 120 km NS and is comprised of 500 sites, which were equipped with different seismic instruments: 3-component 4.5 Hz geophones, short-period Mark 3D instruments, and broadband seismometers (Trillium Compact 120). Details about this deployment are given by Dahm and The Eifel Large-N team (2023).
<table><caption>**Table 1** Search parameters and detection performance for the three study areas in Reykjavíks Peninsula, Iceland; Eifel Volcanic Region, Germany; and Utah FORGE, USA. The searches were conducted on a single compute node with 48 CPU cores, 128 GB RAM, and an Nvidia A100, A40, or A30 GPU for ML inference.</caption><thead><tr><th></th><th>RP, Iceland</th><th>EVR, Germany</th><th>Utah FORGE, USA</th></tr></thead><tbody><tr><td><strong>Dataset</strong></td><td>RP, Iceland</td><td>EVR, Germany</td><td>Utah FORGE, USA</td></tr><tr><td>Seismic stations</td><td>28</td><td>500</td><td>5</td></tr><tr><td>Dataset coverage</td><td>9 months</td><td>1 year</td><td>25 days</td></tr><tr><td>Dataset size</td><td>1 TB</td><td>5 TB</td><td>62 GB</td></tr><tr><td>Sequential waveform duration</td><td>63 years</td><td>1000 years</td><td>125 days</td></tr><tr><td>Processing duration</td><td>17 hours</td><td>2 days</td><td>2.5 hours</td></tr><tr><td><strong>Search Parameters</strong></td><td></td><td></td><td></td></tr><tr><td>ML network; model</td><td>PhaseNet; original</td><td>PhaseNet; original</td><td>PhaseNet; original</td></tr><tr><td>Search area</td><td>62 x 32 km</td><td>60 x 60 km</td><td>9.6 x 6 km</td></tr><tr><td>Search depth</td><td>12 km</td><td>50 km</td><td>4 km</td></tr><tr><td>Search volume</td><td>24 000 km³</td><td>180 000 km³</td><td>230.4 km³</td></tr><tr><td>Octree root node size</td><td>2 km</td><td>2 km</td><td>400 m</td></tr><tr><td>Octree levels</td><td>5 (125 m)</td><td>5 (125 m)</td><td>5 (25 m)</td></tr><tr><td>Velocity model</td><td>Layered 1D</td><td>Layered 1D</td><td>3D</td></tr><tr><td>Magnitude model</td><td>M<sub>w</sub> Dahm et al. (2024)</td><td>ML Stange (2006)</td><td>-</td></tr><tr><td>Velocity model</td><td>Hrubcová and Vavryčuk (2023)</td><td>EDSW, Fig. S19</td><td>Finger et al. (2024)</td></tr><tr><td>Input stretching S<sub>i</sub></td><td>-</td><td>-</td><td>4</td></tr><tr><td>Semblance Sampling Rate</td><td>100 Hz</td><td>100 Hz</td><td>200 Hz</td></tr><tr><td>Semblance threshold</td><td>0.2</td><td>0.2</td><td>0.1</td></tr><tr><td>DeepDenoiser</td><td>No</td><td>No</td><td>No</td></tr><tr><td><strong>Results</strong></td><td></td><td></td><td></td></tr><tr><td>No. qseek catalog</td><td>51 465</td><td>2483</td><td>1464</td></tr><tr><td>No. reference catalog</td><td>18 225</td><td>198</td><td>31</td></tr></tbody></table>

**Analysis** The obtained earthquake catalog, which features SSST corrections, is presented in Figure 10, and the event location uncertainties are presented in Figure S18. Catalogs without corrections and simple SST corrections are presented in Figures S15 and S16. As a benchmark reference for our method, we use the official automatically-generated catalog, which is curated by the Geological Survey of Rhineland Palatinate (LGB) / Erdbebendienst Südwest (EDSW), Freiburg, Germany and based on high-quality permanent seismic stations of the state seismic network (network LE; Erdbebendienst Südwest Baden-Württemberg and Rheinland-Pfalz, 2009) and is shown in Figure S14. Maps of SST delay times for the large-N network are presented in supplementary Figure S1 and S2.

## 3.3 Utah FORGE, USA

Utah FORGE (Frontier Observatory for Research in Geothermal Energy) in southwest Utah, USA, is a field-scale laboratory for studying Enhanced Geothermal Systems (EGS) in granitoid basement rock (Moore et al., 2019). EGS enhances reservoir conductivity at depths typically >2 km by creating fluid pathways in tight basement rock through stimulation, such as hydraulic fracturing. Fluids are injected into an injection well, heated in the fractured reservoir, and extracted via a production well. Key phases include a 2022 stimulation phase creating fluid pathways and a 2023 circulation phase testing the well connection (Moore et al., 2023; McLenan et al., 2023). Utah FORGE, located at the edge of a sedimentary basin overlying dipping basement rock,

requires a 3D velocity model to precisely locate the induced microseismicity due to the sharp velocity contrast.

**Analysis** The earthquake detections take advantage of the 3D ray tracing for travel time calculation using the fast-marching Eikonal solver. We analyze four weeks of seismic data during the circulation tests (July 4 to 31, 2023) recorded by 5 (near-)surface seismometers from the University of Utah Seismograph Stations (UUSS; University of Utah, 1962). The induced microseismicity was similar in event rate and magnitude compared to the 2022 reservoir stimulations, with bursts of hundreds of events in a few hours up to magnitude 0.5 (Niemz et al., 2024). We compare our results to the authoritative regional catalog provided by UUSS and to the relative relocations presented by Niemz et al. (2024). The detections within a 230 km³ source volume revealed induced microseismic events within the geothermal reservoir and natural swarm activity in the Mineral Mountains (Petersen and Pankow, 2023). Details on the circulation-induced microseismicity are presented by Niemz et al. (2024). Detection results using SSSTs are shown in Fig. 11. Plain and SST-corrected versions of the catalog are provided in the supplementary Figure S23.

## 4 Discussion and Conclusions

The developed stacking and migration method and software framework successfully combine machine learning phase annotations and efficient stacking to associate seismic phase arrivals. This facilitates detecting and lo-
<figure><img src="image_11.png" /><figcaption>**Figure 9** Locations of earthquakes on the Reykjanes Peninsula, Iceland detected by qseek with applied SSST corrections. Featuring 50k+ seismic events detected between January 1 and November 1 2020 inside a search volume of 26 000 km³. The high number of events cluster densely along the active structures.</figcaption></figure>

calizing seismic events on an adaptive octree grid with
high precision and resolution. The extraction of sta-
tion corrections through SST and SSST methods further
improves the detection and localization performance.
We demonstrate this in the different analyzed datasets,
where we can increase the number of earthquake detec-
tions (Tab. 2 and S2) and improve the hypocenter loca-
tion certainty (Fig. S8 and S9).

In the volcano-tectonic dataset recorded on Reyk-
janes Peninsula, Iceland, we can identify a notable
number of 51k+ seismic events within the search vol-
ume. The IMO catalog had previously documented
about 18k seismic events from a similar network con-
figuration within the confines of the identical volume.
The seismicity images the geometry of a dike structure
Sigmundsson et al. (2024) with higher fidelity and de-
tail compared with the already high-quality IMO cata-
log. The dense spatial clustering of the events illum-
nates the active geological structures during the unrest
(see supplement for a link to the animation), thus allow-
ing for a more detailed and more constrained interpre-
tation of the seismicity. The location uncertainty is low
in areas with complete azimuthal station coverage and
degrades when the illumination is narrow. The location
uncertainties improve with better station coverage and
the application of SSST corrections (Fig. S8 and Tab. S2).
It should be pointed out that the event depths are sys-

tematically deeper and show a horizontal artifact at ~5.2
km depth. This stems from a velocity contrast within
the layered velocity model (Fig. S13). Flóvenz et al.
(2022) detected 39k+ events within a similar volume and
period, more than detected in this study in the compet-
ing work. However, supplementary DAS data was used
to detect seismicity, resulting in more events and false
detections. A detailed geological interpretation of the
observed seismicity and inverted source mechanisms is
given in Büyükakpinar et al. (2024), Flóvenz et al. (2022)
and Sigmundsson et al. (2024).

In the EVR dataset, we detected 780 events, exceed-
ing the 210 events reported by the regional network
(EDSW). However, this study analyzed seismic data
from the dense local large-N temporary seismic net-
work. In contrast, the EDSW catalog relies on the
more sparse, high-quality stations from the permanent
network. The regional catalog shows tectonic activ-
ity on a local fault system (ODFZ) spatially tightly con-
fined in the qseek catalog. The detections also re-
veal two seismic clusters close to Rieden and Ahrbrück,
shown in the two inlets. The presented seismic cat-
alog qualitatively paints a more complete and precise
picture of the region's activity. The accuracy of loca-
tions is also reflected in the detection and localization
of quarry blasts shown in Figure S17 and the horizon-
tal uncertainty estimations in Figure S8. Further, the
<figure><img src="image_12.png" /><figcaption>**Figure 10** Locations of earthquakes detected by qseek in the Eifel volcanic region between September 20, 2022 and September 1, 2023. A total of 780 events were detected, and the depth color scaling is identical to Fig. S14. The two circular inlets focus on the detected swarms in Altenahr and Rieden. The detections cluster along the main seismic features in the area. SST station corrections were applied.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11** The microseismicity detected at Utah FORGE, USA, during the reservoir circulation tests in July 2023 shows the induced activity within the reservoir and natural swarm activity in the Mineral Mountains to the east. Results from qseek use a local 3D velocity model and SSST corrections containing 1464 event detections.</figcaption></figure>

presented method delivers a continuous resolution of the hypocenter depths. The presented local magnitudes are based on the same $M_L$ model (Stange, 2006), and we obtain similar magnitudes for the same set of events. Dahm and The Eifel Large-N team (2023) presents a geological interpretation of the results and the located DLF events of the large-N network.

Examining seismic data from Utah FORGE, we show that the neural network can annotate induced micro-

seismicity down to $M < -1$ (Fig. S22). This is accomplished by stretching the input waveforms $D(t)$ for the NN and improving the network's annotation sensitivity for microseismic events. Further, the high semblance sampling rate of 200 Hz allows the refinement of the octree down to nodes with 25 m edge length and enhances the precise localization of the seismicity. When comparing the absolute locations from qseek with manually picked relative locations obtained through GlowCrust (Trugman et al., 2022; Niemz et al., 2024), we can chart the lobes and extent of the seismicity cloud around the well strains (Fig. S20). Magnitudes are not presented for the Utah FORGE dataset, but the extraction of spectral magnitudes is planned for future work.

Applying station correction terms in the form of SST and SSST improves the semblance and localization accuracy (Tab. 2 and S2). Incorporating SST and SSST station correction terms refines and further enhances the catalogs by clustering seismic activity onto active structures. These corrections are essential, particularly when only layered 1D velocity models are available. We can demonstrate this effect in our SSST and SST corrected localization when compared to the plain detections presented in the supplementary material. The well-behaved statistical distribution of delay times (Fig. 6) further supports that the bias of travel times is often systematic for individual stations and influenced by the shallow subsurface (Lin and Shearer, 2005). A comparison of the time delay weighting (eq. 10) reveals that the weighted average results in a more significant number
of detected events and a larger cumulative semblance. Conversely, median weighted residuals yield a smaller number of events but at a higher maximum semblance (Tab. S1).

The achieved computational performance of the developed framework surpasses existing open-source alternatives for stacking and migration event detection by a significant margin, and the processing time is reduced from weeks (Lassie) to hours (qseek, Tab. 1). The performance enhancement is rooted in the octree search and resolute application of asynchronous programming principles throughout the advanced software stack. The computational speed-up is critical for processing large seismic datasets and iterative computation of seismic catalogs. It enables the generation of robust, complete, and precise seismic catalogs from large networks and across extensive regions. Utilizing an NN image function instead of the STA/LTA image functions stabilizes and simplifies the stacking and migration approach: it eliminates the need to normalise STA/LTA counts (raw waveforms), resulting in a more robust detection of events. The shape of the NN annotation function propagates phase and arrival uncertainties into the backprojection stack. This, together with the adaptive refinement of the octree, leads to more precise hypocenter locations. However, the fine resolution of the octree grid can lead to high memory consumption and require capable computing hardware. This also makes the search within large volumes inefficient. These shortcomings can be counteracted by implementing the stacking and migration algorithm (parstack) in single precision floats (16-bit), effectively cutting the memory consumption in half and doubling the processing speed. This is planned for future work. We encourage contributing to the open-source software framework. To foster scientific collaboration and ease the development experience we designed the framework in a modular fashion and adhere to modern open collaboration platforms (i.e. GitHub), code formatting, and linting, which eases on-boarding of scientists and developers alike.

Our approach of combining machine learning with stacking and migration on an adaptive octree grid has identified a large number of seismic events with high spatial resolution in different geological settings and seismic datasets. Accurate and complete catalogs can offer a richer understanding of the geological and seismo-tectonic dynamics beyond what standard location catalogs may provide. By creating qseek, we want to enable the analysis of large datasets. This facilitates new perspectives and insights for discoveries in Earth sciences.

## Acknowledgements

We thank three anonymous reviewers for taking the time and effort to review the manuscript. We sincerely appreciate all the valuable comments and suggestions that helped us improve the manuscript's quality. We thank and acknowledge the recording and sharing of seismic data by the Iceland Meteorological Office, the University of Iceland, and colleagues from the Czech Academy of Sciences, as well as colleagues from GFZ

<table><caption>**Table 2** Event detection and semblance (E) statistics with and without station corrections. Station corrections lead to more event detections and increased cumulative and maximum semblance.</caption><thead><tr><th>Corrections</th><th>None</th><th>SST</th><th>SSST</th></tr></thead><tbody><tr><td><strong>RP, Iceland</strong></td><td></td><td></td><td></td></tr><tr><td>Detections</td><td>50 744</td><td>50 826</td><td>51 465</td></tr><tr><td>E Cum.</td><td>28 193</td><td>29 794</td><td>31 725</td></tr><tr><td>E Avg.</td><td>0.555</td><td>0.586</td><td>0.616</td></tr><tr><td>E Max.</td><td>1.499</td><td>1.641</td><td>1.774</td></tr><tr><td><strong>EVR, Germany</strong></td><td></td><td></td><td></td></tr><tr><td>Detections</td><td>2362</td><td>2364</td><td>2483</td></tr><tr><td>E Cum.</td><td>1 177.1</td><td>1 319.6</td><td>1 483.1</td></tr><tr><td>E Avg.</td><td>0.259</td><td>0.297</td><td>0.312</td></tr><tr><td>E Max.</td><td>0.958</td><td>1.190</td><td>1.351</td></tr><tr><td><strong>Utah FORGE, USA</strong></td><td></td><td></td><td></td></tr><tr><td>Detections</td><td>1431</td><td>1448</td><td>1464</td></tr><tr><td>E Cum.</td><td>602.062</td><td>610.805</td><td>615.781</td></tr><tr><td>E Avg.</td><td>0.421</td><td>0.422</td><td>0.421</td></tr><tr><td>E Max.</td><td>1.443</td><td>1.400</td><td>1.420</td></tr></tbody></table>

Claus Milkereit and Stefan Mikulla installing the MAGIC seismic network during the 2020 unrest on the Reykjavanes Peninsula. Further, we want to express our acknowledgment to the GFZ Helmholtz-Zentrum for funding the GIPP Pool and the Eifel Large-N deployment, as well as supporting the MAGIC HART activities and enabling the swift scientific reaction to the seismic crisis in Iceland. P.N. was funded by the U.S. DOE under grant DE-EE0007080 "Enhanced Geothermal System Concept Testing and Development at the Milford City, Utah FORGE Site". P. B. is funded by project 517982028 (BU 4346/1-1) of the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation). This work utilized high-performance computing resources made possible by funding from the Ministry of Science, Research and Culture of the State of Brandenburg (MWFK). It is operated by the IT Services and Operations unit of the Helmholtz Centre Potsdam.

## Data and code availability

The software framework qseek is available as open-source software at https://github.com/pyrocko/qseek. Active participation in code reviews, software development, and framework extension is encouraged. Online documentation and set-up instructions are available at https://pyrocko.github.io/qseek/.

Analysis and plots are generated using the framework and scripts using Matplotlib (Hunter, 2007), and maps were created using QGIS. Earthquake catalogs and velocity models are available on Zenodo (Isken, 2025). Seismic data recorded by MAGIC Project (9H; Dahm et al., 2020a) is available through the GEOFON FDSN WebServices available at https://geofon.gfzpotsdam.de/waveform/webservices/fdsnws.php. Other restricted data was recorded by the Reykjanet (network code 7E; Josef Horalek, 2013) and the Iceland National Network (SI; Icelandic Meteorological Office, 1992). The IMO catalog for Reykjanes Peninsula was
provided by personal communication. The EDSW cat-
alog for the Eifel Volcanic Region is accessible through
WFS services at https://mapserver.lgb-rlp.de/cgi-bin/mc_
erdbeben. Utah FORGE relative relocations are available
from Niemz et al. (2024).

Competing interests

The authors declare that they have no conflict of inter-
est.

References

Allen, R. Automatic phase pickers: Their present use and future prospects. *Bulletin of the Seismological Society of America*, **72** (6B):S225-S242, 1982.

Bagagli, M. filter picker: general purpose, broad-band, phase de-
tector and picker, May 2022. doi: 10.5281/zenodo.6511985.

Beaucé, E., Frank, W. B., Seydoux, L., Poli, P., Groebner, N., van der Hilst, R. D., and Campillo, M. BPMF: A Backprojection and Matched-Filtering Workflow for Automated Earthquake Detection and Location. *Seismological Research Letters*, 95(2A): 1030–1042, 12 2023. doi: 10.1785/0220230230.

Billings, S. D., Sambridge, M. S., and Kennett, B. L. N. Errors in hypocenter location: Picking, model, and magnitude dependence. *Bulletin of the Seismological Society of America*, 84(6): 1978–1990, 12 1994. doi: 10.1785/BSSA0840061978.

Bormann, P. and Dewey, J. W. The new IASPEI standards for determining magnitudes from digital data and their relation to classical magnitudes. In *New manual of seismological observatory practice 2 (NMSOP-2)*, pages 1–44. Deutsches Geo-ForschungsZentrum GFZ, 2012. doi: 10.2312/GFZ.NMSOP-2_IS_3.3.

Büyükakpınar, P., Isken, M. P., Heimann, S., Dahm, T., Kühn, D., Starke, J., López Comino, J. A., Cesca, S., Doubravová, J., Gudnason, E. A., and Ágústsdóttir, T. Understanding the Seismic Signature of Transtensional Opening in the Reykjanes Peninsula Rift Zone, SW Iceland. *Journal of Geophysical Research: Solid Earth*, 2024. doi: 10.1029/2024jb029566.

Cesca, S. and Grigoli, F. Chapter Two - Full Waveform Seismological Advances for Microseismic Monitoring, volume 56 of Advances in Geophysics. Elsevier, 2015. doi: 10.1016/bs.agph.2014.12.002.

Cesca, S., Letort, J., Razafindrakoto, H. N. T., Heimann, S., Rivalta, E., Isken, M. P., Nikkhoo, M., Passarelli, L., Petersen, G. M., Coton, F., and Dahm, T. Drainage of a deep magma reservoir near Mayotte inferred from seismicity and deformation. *Nature Geoscience*, 13(1):87–93, Jan. 2020. doi: 10.1038/s41561-019-0505-5.

Dahm, T. and The Eifel Large-N team. A large-N passive seismological experiment to unravel the structure and activity of the transtrustal magma system of the Eifel Volcanic Field. *EGU General Assembly, Vienna, Austria*, 2023. doi: 10.5194/egusphere-egu23-2590.

Dahm, T., Jousset, P., Heimann, S., Milkereit, C., Hersir, G. P., and Magnússon, R. MAGIC - Seismic network MAGma in ICeland, 2020a. doi: 10.14470/4U7575229166.

Dahm, T., Stiller, M., Mechie, J., Heimann, S., Hensch, M., Woith, H., Schmidt, B., Gabriel, G., and Weber, M. Seismological and Geophysical Signatures of the Deep Crustal Magma Systems of the Cenozoic Volcanic Fields Beneath the Eifel, Germany. *Geochemistry, Geophysics, Geosystems*, 21(9), 2020b. doi: 10.1029/2020GC009062.

Dahm, T., Kühn, D., Cesca, S., Isken, M., and Heimann, S. Earth-
quake Moment Magnitudes from Peak Ground Displacements
and Synthetic Green's Functions. *Seismica*, 3(2), Dec. 2024. doi:
10.26443/seismica.v3i2.1205.

Defourny, A., Blard, P.-H., Zimmermann, L., Jobé, P., Collignon, A., Nguyen, F., and Dassargues, A. δ¹³C, CO₂ /³He and ³He / ⁴He ratios reveal the presence of mantle gas in the CO₂-rich groundwatersof the Ardennes massif (Spa, Belgium). *Hydrology and Earth System Sciences*, 26(10):2637–2648, 2022. doi: 10.5194/hess-26-2637-2022.

del Fresno, C., Cesca, S., Klügel, A., Domínguez Cerdeña, I., Díaz-Suárez, E. A., Dahm, T., García-Cañada, L., Meletlidis, S., Milkereit, C., Valenzuela-Malebrán, C., López-Díaz, R., and López, C. Magmatic plumbing and dynamic evolution of the 2021 La Palma eruption. *Nature Communications*, 14(1), Jan. 2023. doi: 10.1038/s41467-023-35953-y.

Dimri, V., Srivastava, R., and Vedanti, N. Chapter 8 - Hamilton's ray equations. In Slawinski, M., editor, *Seismic Waves and Rays in Elastic Media*, volume 34 of *Handbook of Geophysical Exploration: Seismic Exploration*, pages 173–195. Pergamon, 2003. doi: 10.1016/S0950-1401(03)80057-5.

Dougherty, S. L., Cochran, E. S., and Harrington, R. M. The LArge-n Seismic Survey in Oklahoma (LASSO) Experiment. *Seismological Research Letters*, 90(5):2051–2057, 08 2019. doi: 10.1785/0220190094.

Drew, J., White, R. S., Tilmann, F., and Tarasewicz, J. Coalescence microseismic mapping. *Geophysical Journal International*, 195 (3):1773–1785, 09 2013. doi: 10.1093/gji/ggt331.

Ducrocq, C., Árnadóttir, T., Einarsson, P., Jónsson, S., Drouin, V., Geirsson, H., and Hjartardóttir, Á. R. Widespread fracture movements during a volcano- tectonic unrest: the Reykjanes Peninsula, Iceland, from 2019–2021 TerraSAR-X interferometry. *Bull. Volcanol*, 86(2), Jan. 2024. doi: 10.1007/s00445-023-01699-0.

Ellsworth, W. L. Injection-Induced Earthquakes. *Science*, 341 (6142):1225942, 2013. doi: 10.1126/science.1225942.

Erdbebendienst Südwest Baden-Württemberg and Rheinland-
Pfalz. Erdbebendienst Südwest, 2009. doi: 10.7914/SN/LE.

Finger, C., Niemz, P., Ermert, L., and Lanza, F. A Composite 3D Velocity Model for Utah FORGE. Technical report, University of Utah Seismograph Stations, 2024. doi: 10.15121/2305384.

Flóvenz, Ó. G., Wang, R., Hersir, G. P., Dahm, T., Hainzl, S., Vasileva, M., Drouin, V., Heimann, S., Isken, M. P., Gudnason, E. Á., et al. Cyclical geothermal unrest as a precursor to Iceland's 2021 Fagradalsfjall eruption. *Nature Geoscience*, 15(5):397–404, 2022. doi: http://dx.doi.org/10.1038/s41561-022-00930-5.

Gal, F., Michel, B., Gilles, B., Frédéric, J., and Karine, M. CO₂ escapes in the Laacher See region, East Eifel, Germany: Application of natural analogue onshore and offshore geochemical monitoring. *International Journal of Greenhouse Gas Control*, 5 (4):1099–1118, 2011. doi: 10.1016/j.ijggc.2011.04.004.

Gharti, H. N., Oye, V., Roth, M., and Kühn, D. Automated mi-
croearthquake location using envelope stacking and robust
global optimization. *Geophysics*, 75(4):MA27-MA46, 08 2010.
doi: 10.1190/1.3432784.

Greenfield, T., White, R. S., Winder, T., and Ágústsdóttir, T. Seismicity of the Askja and Bárðarbunga volcanic systems of Iceland, 2009–2015. *Journal of Volcanology and Geothermal Research*, 391:106432, 2020. doi: 10.1016/j.jvolgeores.2018.08.010.

Greenfield, T., Winder, T., Rawlinson, N., Maclennan, J., White, R. S., Ágústsdóttir, T., Bacon, C. A., Brandsdóttir, B., Eibl, E. P. S., Glastonbury-Southern, E., Gudnason, E. A., Hersir, G. P., and Horálek, J. Deep long period seismicity preceding and during the 2021 Fagradalsfjall eruption, Iceland. *Bulletin of Volcanology*, 84(12), Nov. 2022. doi: 10.1007/s00445-022-01603-2.
Grigoli, F., Cesca, S., Vassallo, M., and Dahm, T. Automated seismic event location by travel-time stacking: An application to mining induced seismicity. *Seismological Research Letters*, 84(4): 666-677, 2013. doi: 10.1785/0220120191.

Grigoli, F., Cesca, S., Amoroso, O., Emolo, A., Zollo, A., and Dahm, T. Automated seismic event location by waveform coherence analysis. *Geophysical Journal International*, 196(3):1742-1753, 2014. doi: 10.1785/0220120191.

Grigoli, F., Cesca, S., Krieger, L., Kriegerowski, M., Gammaldi, S., Horalek, J., Priolo, E., and Dahm, T. Automated microseismic event location using master-event waveform stacking. *Scientific reports*, 6(1):1-13, 2016. doi: 10.1038/srep25744.

Hamilton, W. R. Theory of systems of rays. *The Transactions of the Royal Irish Academy*, pages 69-174, 1828.

Heimann, S., Kriegerowski, M., Isken, M., Cesca, S., Daout, S., Grigoli, F., Juretzek, C., Megies, T., Nooshiri, N., Steinberg, A., et al. *Pyrocko-An open-source seismology toolbox and library*, 2017. doi: 10.5880/GFZ.2.1.2017.001.

Heimann, S., Vasyura-Bathke, H., Sudhaus, H., Isken, M. P., Kriegerowski, M., Steinberg, A., and Dahm, T. A Python framework for efficient use of pre-computed Green's functions in seismological and other physical forward and inverse source problems. *Solid Earth*, 10(6):1921-1935, 2019. doi: 10.5194/se-10-1921-2019.

Hensch, M., Dahm, T., Ritter, J., Heimann, S., Schmidt, B., Stange, S., and Lehmann, K. Deep low-frequency earthquakes reveal ongoing magmatic recharge beneath Laacher See Volcano (Eifel, Germany). *Geophysical Journal International*, 216(3): 2025-2036, 01 2019. doi: 10.1093/gji/ggy532.

Hinzen, K.-G. Stress field in the Northern Rhine area, Central Europe, from earthquake fault plane solutions. *Tectonophysics*, 377(3):325-356, 2003. doi: 10.1016/j.tecto.2003.10.004.

Horálek, J., Fischer, T., Boušková, A., and Jedlička, P. Western Bohemia/Vogtland in the light of the WEBNET network. *Stud. Geophys. Geod.*, 44(2):107-125, 2000. doi: 10.1023/A:1022198406514.

Hrubcová, P. and Vavryčuk, V. Tectonic stress changes related to plate spreading prior to the 2021 Fagradalsfjall eruption in SW Iceland. *Tectonophysics*, 851:229761, 2023. doi: https://doi.org/10.1016/j.tecto.2023.229761.

Hunter, J. D. Matplotlib: A 2D graphics environment. *Computing in Science & Engineering*, 9(3):90-95, 2007. doi: 10.1109/M-CSE.2007.55.

Hutton, L. K. and Boore, D. M. The ML scale in Southern California. *Bulletin of the Seismological Society of America*, 77(6): 2074-2094, Dez 1987.

Icelandic Meteorological Office. Icelandic National Digital Seismograph Network, 1992. doi: 10.7914/ZKJM-TJ71.

Isken, M. Supplement to Qseek: A Data-Driven Framework for Automated Earthquake Detection, Localization and Characterization, 2025. doi: 10.5281/zenodo.14802237.

Josef Horalek. Reykjanet, 2013. doi: 10.7914/SN/7E_2013.

Kao, H. and Shan, S.-J. The source-scanning algorithm: Mapping the distribution of seismic sources in time and space. *Geophysical Journal International*, 157(2):589-594, 2004. doi: 10.1111/j.1365-246X.2004.02276.x.

Li, L., Tan, J., Schwarz, B., Staněk, F., Poiata, N., Shi, P., Diekmann, L., Eisner, L., and Gajewski, D. Recent Advances and Challenges of Waveform-Based Seismic Location Methods at Multiple Scales. *Reviews of Geophysics*, 58(1):e2019RG000667, 2020. doi: 10.1029/2019RG000667.

Lin, G. and Shearer, P. Tests of relative earthquake location techniques using synthetic data. *Journal of Geophys-

ical Research: Solid Earth, 110(B4), 2005. doi: https://doi.org/10.1029/2004JB003380.

Lomax, A. and Savvaidis, A. High-Precision Earthquake Location Using Source-Specific Station Terms and Inter-Event Waveform Similarity. *Journal of Geophysical Research: Solid Earth*, 127(1): e2021JB023190, 2022. doi: 10.1029/2021JB023190.

Lomax, A., Virieux, J., Volant, P., and Berge-Thierry, C. *Probabilistic Earthquake Location in 3D and Layered Models*, page 101-134. Springer Netherlands, 2000. doi: 10.1007/978-94-015-9536.

Lomax, A., Satriano, C., and Vassallo, M. Automatic Picker Developments and Optimization: FilterPicker—a Robust, Broadband Picker for Real-Time Seismic Monitoring and Earthquake Early Warning. *Seismological Research Letters*, 83(3):531-540, 05 2012. doi: 10.1785/gssrl.83.3.531.

McLennan, J., England, K., Rose, P., Moore, J., and Barker, B. Stimulation of a High-Temperature Granitic Reservoir at the Utah FORGE Site. In *Day 3 Thu, February 02, 2023*, page D031S008R001. SPE, 2023. doi: 10.2118/212346-MS.

Meagher, D. Geometric modeling using octree encoding. *Computer Graphics and Image Processing*, 19(2):129-147, 1982. doi: 10.1016/0146-664X(82)90104-6.

Mertes, H. *Aufbau und Genese des Westeifeler Vulkanfeldes*. Number 9 in 315. Institut für Geologie der Ruhr-Universität-Bochum, 1983.

Moore, J., McLENNAN, J., Allis, R., Pankow, K., Simmons, S., Podgorney, R., Wannamaker, P., and Rickard, W. The Utah Frontier Observatory for Research in Geothermal Energy (FORGE): An International Laboratory for Enhanced Geothermal System Technology Development. In *PROCEEDINGS*, pages SGP-TR-214, Stanford University, Stanford, California, USA, 2019.

Moore, J., McLennan, J., Pankow, K., Finnila, A., Dyer, B., Karvounis, D., Bethmann, F., Podgorney, R., Rutledge, J., Meir, P., Xing, P., Jones, C., Barker, B., Simmons, S., and Damjanac, B. Current Activities at the Utah Frontier Observatory for Research in Geothermal Energy (FORGE): A Laboratory for Characterizing, Creating and Sustaining Enhanced Geothermal Systems. In *All Days*, pages ARMA-2023-0749. ARMA, 2023. doi: 10.56952/ARMA-2023-0749.

Mousavi, S. M., Ellsworth, W. L., Zhu, W., Chuang, L. Y., and Beroza, G. C. Earthquake transformer—an attentive deep-learning model for simultaneous earthquake detection and phase picking. *Nature Communications*, 11(1):1-12, 2020. doi: 10.1038/s41467-020-17591-w.

Münchmeyer, J. PyOcto: A high-throughput seismic phase associator. *Seismica*, 3(1), Jan. 2024. doi: 10.26443/seismica.v3i1.1130.

Münchmeyer, J., Woollam, J., Rietbrock, A., Tilmann, F., Lange, D., Bornstein, T., Diehl, T., Giunchi, C., Haslinger, F., Jozinović, D., Michelini, A., Saul, J., and Soto, H. Which Picker Fits My Data? A Quantitative Evaluation of Deep Learning Based Seismic Pickers. *Journal of Geophysical Research: Solid Earth*, 127 (1):e2021JB023499, 2022. doi: 10.1029/2021JB023499.

Niemz, P., McLennan, J., Pankow, K. L., Rutledge, J., and England, K. Circulation experiments at Utah FORGE: Near-surface seismic monitoring reveals fracture growth after shut-in. *Geothermics*, 119:102947, 2024. doi: 10.1016/j.geothermics.2024.102947.

Nooshiri, N., Saul, J., Heimann, S., Tilmann, F., and Dahm, T. Revision of earthquake hypocentre locations in global bulletin data sets using source-specific station terms. *Geophysical Journal International*, 208(2):589-602, 10 2016. doi: 10.1093/gji/ggw405.

Obermann, A., Sánchez-Pastor, P., Wu, S., Wollin, C., Baird, A. F., Isken, M. P., Clinton, J., Goertz-Allmann, B. P., Dahm, T., Wuestefeld, A., Shi, P., Lanza, F., Gyger, L., Wetter, S., Hjörleifsdóttir, V., Langet, N., Brynjarsson, B., Jousset, P., and Wiemer, S. Com-
bined Large-N Seismic Arrays and DAS Fiber Optic Cables across the Hengill Geothermal Field, Iceland. *Seismological Research Letters*, 93(5):2498-2514, 07 2022. doi: 10.1785/0220220073.

Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. Pytorch: An imperative style, high-performance deep learning library. *Advances in neural information processing systems*, 32, 2019.

Petersen, G. M. and Pankow, K. L. Small-Magnitude Seismic Swarms in Central Utah (US): Interactions of Regional Tectonics, Local Structures and Hydrothermal Systems. *Geochemistry, Geophysics, Geosystems*, 24(7):e2023GC010867, 2023. doi: 10.1029/2023GC010867.

Poiata, N., Satriano, C., Vilotte, J.-P., Bernard, P., and Obara, K. Multiband array detection and location of seismic sources recorded by dense seismic networks. *Geophysical Journal International*, 205(3):1548-1573, 02 2016. doi: 10.1093/gji/ggw071.

Reinig, F., Wacker, L., Jöris, O., Oppenheimer, C., Guidobaldi, G., Nievergelt, D., Adolphi, F., Cherubini, P., Engels, S., Esper, J., et al. Precise date for the Laacher See eruption synchronizes the Younger Dryas. *Nature*, 595(7865):66-69, 2021. doi: 10.1038/s41586-021-03608-x.

Richards-Dinger, K. B. and Shearer, P. M. Earthquake locations in southern California obtained using source-specific station terms. *Journal of Geophysical Research: Solid Earth*, 105(B5):10939-10960, 2000. doi: https://doi.org/10.1029/2000JB900014.

Ross, Z. E., Meier, M., Hauksson, E., and Heaton, T. H. Generalized Seismic Phase Detection with Deep Learning. *Bulletin of the Seismological Society of America*, 108(5A):2894-2901, 08 2018. doi: 10.1785/0120180080.

Sethian, J. A. A fast marching level set method for monotonically advancing fronts. *Proceedings of the National Academy of Sciences*, 93(4):1591-1595, 1996. doi: 10.1073/pnas.93.4.1591.

Shi, P., Grigoli, F., Lanza, F., Beroza, G. C., Scarabello, L., and Wiemer, S. MALMI: An Automated Earthquake Detection and Location Workflow Based on Machine Learning and Waveform Migration. *Seismological Research Letters*, 93(5):2467-2483, 05 2022. doi: 10.1785/0220220071.

Sigmundsson, F., Parks, M., Geirsson, H., Hooper, A., Drouin, V., Vogfjörð, K. S., Ófeigsson, B. G., Greiner, S. H. M., Yang, Y., Lanzi, C., Pascale, G. P. D., Jónsdóttir, K., Hreinsdóttir, S., Tolpekin, V., Friðriksdóttir, H. M., Einarsson, P., and Barsotti, S. Fracturing and tectonic stress drives ultrarapid magma flow into dikes. *Science*, 1(1):eadn2838, 2024. doi: 10.1126/science.adn2838.

Stange, S. ML Determination for Local and Regional Events Using a Sparse Network in Southwestern Germany. *Journal of Seismology*, 10(2):247-257, Apr. 2006. doi: 10.1007/s10950-006-9010-6.

Stein, S. and Wysession, M. *An introduction to seismology, earthquakes, and earth structure*. Blackwell Science, Philadelphia, PA, Sept. 2002.

Trugman, D. T., Chamberlain, C. J., Savvaidis, A., and Lomax, A. GrowClust3D.Jl: A Julia Package for the Relative Relocation of Earthquake Hypocenters Using 3D Velocity Models. *Seismological Research Letters*, 94(1):443-456, 2022. doi: 10.1785/0220220193.

Uhrhammer, R. A. and Collins, E. R. Synthesis of Wood-Anderson seismograms from broadband digital records. *Bulletin of the Seismological Society of America*, 80(3):702-716, 06 1990. doi: 10.1785/BSSA0800030702.

University of Utah. University of Utah Regional Seismic Network, 1962. doi: 10.7914/SN/UU.

Vasyura-Bathke, H., Dettmer, J., Biegel, K., Salvage, R. O., Eaton, D., Ackerley, N., Samsonov, S., and Dahm, T. Bayesian inference

elucidates fault-system anatomy and resurgent earthquakes induced by continuing saltwater disposal. *Communications Earth and Environment*, 4(1), Nov. 2023. doi: 10.1038/s43247-023-01064-1.

Wilding, J. D., Zhu, W., Ross, Z. E., and Jackson, J. M. The magmatic web beneath Hawai'i. *Science*, 379(6631):462-468, 2023. doi: 10.1126/science.ade5755.

Winder, T., Bacon, C., Smith, J. D., Hudson, T. S., Drew, J., and White, R. S. *QuakeMigrate v1.0.0*, Jan. 2021. doi: 10.5281/zenodo.4442749.

Woollam, J., Münchmeyer, J., Tilmann, F., Rietbrock, A., Lange, D., Bornstein, T., Diehl, T., Giunchi, C., Haslinger, F., Jozinović, D., et al. SeisBench—A toolbox for machine learning in seismology. *Seismological Research Letters*, 93(3):1695-1709, 2022.

Zhu, W. and Beroza, G. C. PhaseNet: a deep-neural-network-based seismic arrival-time picking method. *Geophysical Journal International*, 216(1):261-273, 10 2018. doi: 10.1093/gji/ggy423.

Zhu, W., Mousavi, S. M., and Beroza, G. C. Seismic signal denoising and decomposition using deep neural networks. *IEEE Transactions on Geoscience and Remote Sensing*, 57(11):9476-9488, 2019. doi: 10.1109/TGRS.2019.2926772.

Zhu, W., McBrearty, I. W., Mousavi, S. M., Ellsworth, W. L., and Beroza, G. C. Earthquake Phase Association Using a Bayesian Gaussian Mixture Model. *Journal of Geophysical Research: Solid Earth*, 127(5), 2022. doi: 10.1029/2021JB023249.

Zolitschka, B., Negendank, J. F., and Lottermoser, B. Sedimentological proof and dating of the early Holocene volcanic eruption of Ulmener Maar (Vulkaneifel, Germany). *Geologische Rundschau*, 84:213-219, 1995.

The article Qseek: A Data-Driven Framework for Automated Earthquake Detection, Localization and Characterization © 2025 by Marius Paul Isken is licensed under CC BY 4.0.