Original Article

Urban challenges in seismology: seismic monitoring of Kwintsheul’s geothermal operation (the Netherlands)

David Naranjo¹, Marius Isken², Boris Boullenger³, Tania Toledo⁴, Cornelis Weemstra^5,1 and Deyan Draganov¹

¹Department of Geoscience and Engineering, Delft University of Technology, Delft, the Netherlands; ²GFZ Helmholtz Centre for Geosciences, Potsdam, Germany; ³TNO Geological Survey of the Netherlands, Utrecht, the Netherlands; ⁴Swiss Seismological Service (SED) at ETH Zurich, Switzerland; ⁵Department of Seismology and Acoustics, Royal Netherlands Meteorological Institute, De Bilt, the Netherlands

Abstract

Seismic monitoring is essential for understanding subsurface processes, particularly in geothermal operations where low-magnitude events can provide valuable insights into reservoir behaviour. There are two significant challenges when monitoring the seismicity in Dutch geothermal operations: (1) detecting signals from seismic events as noise levels are typically high in regions hosting geothermal operations, and (2) accurately estimating their corresponding hypocentre and uncertainty. In this study, we present a comprehensive workflow for detecting and characterising low-magnitude seismic events. Specifically, we integrated data preparation, template-matching and machine-learning-based event detection, and probabilistic hypocentre estimation. Applying this workflow to 4 months of recordings in Kwintsheul, Netherlands, we detected 65 events with coherent signals, including six weak seismic events (ML < 0.0) near a local fault and a geothermal injection well. These events suggest the presence of a recurring microseismic sequence previously unreported in the area. However, spatial uncertainties, the short monitoring period, and the limited azimuthal coverage make the nature of these events unclear. Our findings highlight the importance of improving network design and refining velocity models to reduce uncertainties in event locations and magnitudes. The proposed workflow offers a scalable solution for enhancing seismic monitoring, particularly in urban and geothermal settings.

Keywords: Seismic-event detection; hypocentre inversion; geothermal energy

 

 

1. Introduction

The Netherlands is increasingly adopting geothermal energy to reduce carbon emissions and transition to local, sustainable energy sources. Understanding the subsurface processes associated with geothermal operations is essential to optimise resource extraction and ensure long-term sustainability. These processes, such as fluid migration and fault reactivation, are often accompanied by low-magnitude seismic events. Therefore, seismic monitoring is a key tool for characterising geothermal reservoirs.

In 2018, a geothermal doublet became operational in Kwintsheul, South Holland, to heat 64 hectares of greenhouses. At that time, the coverage of the Royal Netherlands Meteorological Institute’s (KNMI) seismological network was still relatively sparse in the province of South Holland, as shown in Figure 1. In 2019, a temporary passive seismic network was installed to study the interaction between the geothermal operation and the underlying geological structures. During this monitoring period, a seismic event with a magnitude of M_d 0.16 was detected on 14 July 2019 (Muntendam-Bos et al., 2022). This event suggested the presence of more low-magnitude seismicity in the region, motivating us to investigate the feasibility of detecting additional events in the recorded data.

Figure 1. (A) Regional distribution of seismicity in the Netherlands, main structural features, and location of the Kwintsheul geothermal doublet. The epicentres from the complete KNMI earthquake catalogue (1920–2021) (KNMI, 2023) are shown in red (natural seismicity) and blue (induced seismicity) circles. The stations of KNMI’s permanent seismic network (operational from 22 July 2019 to 9 November 2019) are shown in green (KNMI, 1993). The West Netherlands Basin (WNB), Groningen region, and Roer Valley Graben (RVG) are highlighted, and the individual faults (v. Gessel et al., 2021) are shown in grey. The orange polygon delineates the province of South Holland, where the Kwintsheul geothermal doublet is located (black square). (B) Enlarged view of the Kwintsheul geothermal site. The injector well (KW-GT-01) is shown in blue, and the producer well (KW-GT-02) is shown in red. The temporal monitoring seismic stations are shown in dark red triangles. These were operational from 22 July to 9 November 2019. The background map tiles are provided by OpenStreetMap (OpenStreetMap Contributors, 2017).

Monitoring low-magnitude seismicity around geothermal operations is challenging in the Netherlands. Many of these operations are located in noisy urban environments, where low-magnitude events are difficult to detect due to strong background noise (Groos & Ritter, 2009). This raises a fundamental question: how can we effectively monitor low-magnitude seismicity in such settings to better characterise subsurface processes?

We address this question by presenting a comprehensive workflow for detecting and characterising low-magnitude seismic events in urban areas. The workflow integrates data management, machine-learning-based event detection, and probabilistic hypocentre estimation, providing a portable approach to seismic monitoring. We apply this workflow to a 4-month monitoring period in Kwintsheul and demonstrate its effectiveness in detecting previously unreported seismicity. Our findings highlight the importance of dense seismic arrays to enhance redundancy and azimuthal coverage, which is crucial for detecting and locating low-magnitude seismicity. The results provide insights into ongoing low-magnitude (M_L < 0) seismic activity; however, the nature of these events remains inconclusive due to limitations in the velocity model and the short monitoring period.

2. Geological context and seismicity in the Netherlands

The complete earthquake catalogue for the Netherlands and surrounding regions (KNMI, 2023) indicates that natural seismicity is concentrated in the Roer Valley Graben (RVG), depicted as red circles in Figure 1. These events are associated with the northeastern extension of the active faults bounding the Lower Rhine Graben (Houtgast & van Balen, 2000). Induced seismic activity is concentrated in the northern Netherlands, particularly in the Groningen region, linked to gas extraction (Muntendam-Bos et al., 2022). These induced events, shown as blue circles in Figure 1, are typically low in magnitude (M_L ≤ 4.0).

Kwintsheul is located on top of the West Netherlands Basin (WNB) in the province of South Holland. To this date, no seismic activity has been detected in this region by KNMI’s seismic monitoring network (KNMI, 2023). The WNB is a 60-km-wide transtensional basin that forms part of a failed rift system (Boersma et al., 2021). It contains NW-SE oriented normal faults and consists of Permian to Tertiary deposits that reach thicknesses of up to 5 km, with a well-connected fault system throughout (Worum et al., 2005; Duin et al., 2006; Boersma et al., 2021). The location of the geothermal doublet is depicted in Figure 1B. The injection and production wells reach a depth of approximately 2300 m, targeting the Delft Sandstone Member. The separation between the two wells is approximately 1500 m at reservoir depth. Until 2020, KNMI’s magnitude of completeness around Kwintsheul was M_L1.0, meaning earthquakes below this threshold were likely undetected (Muntendam-Bos et al., 2022; Ruigrok et al., 2023).

3. Workflow for seismic monitoring

To characterise seismicity around Kwintsheul’s geothermal operation, we adopted a multi-stage workflow designed to manage seismic data, detect seismic events, and characterise their hypocentres and magnitudes (Figure 2). The workflow begins with data preparation, including seismic waveform management, station information, and retrieving local velocity models. Next, event detections are performed using a combination of template matching (TM) and machine-learning-based detection techniques. Finally, the detected events are characterised by estimating their hypocentres and magnitudes, accounting for uncertainties in both observed and theoretical phase arrival times. The following sections provide a detailed description of each stage in the workflow.

Figure 2. Workflow for seismic-event detection, location, and characterisation. The process is divided into four main stages: data management, event detection, hypocentre estimation, and local magnitude estimation.

3.1 Data management

Waveform data

Our workflow begins with preparing seismic waveform data and station metadata. We use the Squirrel package from the Pyrocko toolbox (Heimann et al., 2017) to manage these datasets. Its metadata caching allows quick inspection, helping us detect and resolve issues such as gaps or timing errors.

Velocity-model retrieval

As part of the data management stage, we also retrieve the velocity model required for hypocentre inversion. In this study, we use Velmod (3.1), the Dutch regional seismic velocity model (Pluymaekers et al., 2017). Velmod is a 3D velocity model that integrates velocities measured in boreholes (sonic logs and check-shot data) with stacking velocities derived from seismic surveys. It provides a 3D distribution of seismic velocities and includes estimates of their uncertainties. Further details regarding the parameterisation of Velmod 3.1 are provided in Section 3.3. To retrieve the local velocity cube from Velmod, we use PRESEIS (Kraaijpoel, 2025), which extracts the region of interest and outputs the model as a multidimensional array in Xarray format (Hoyer & Hamman, 2017). The next stage of the workflow involves event detection and preliminary locations.

3.2 Event detection

Template-matching detection – EQcorrscan

As the waveforms of a previously detected seismic event are available, we use TM to identify additional events with similar waveform characteristics. In this method, previously detected events serve as templates, which are cross-correlated with continuous waveform recordings from multiple seismic stations to find matching events. Events identified through this process are known as repeating microearthquakes, as they exhibit high waveform similarity across different stations (e.g. Ellsworth, 1995; Menke, 1999). These events often originate from a concentrated volume of hypocentres and provide valuable insights for monitoring induced seismicity.

TM was performed using the open-source Python package EQcorrscan (Chamberlain et al., 2017). The preprocessing parameters and threshold values for the TM are presented in Section 5.1.

Machine-learning event detection

Machine-learning event detection involves two steps. Firstly, we identify the onset times of seismic phases, such as P- and S-waves, by distinguishing their unique signal attributes from background noise. Secondly, we associate these seismic phases across multiple stations with a single seismic event. These steps are explained in detail below and the preprocessing parameters are introduced in Section 5.

Step 1: Waveform image function

Deep-learning phase-picking models are trained on large datasets of manually labelled seismic phases (Mousavi et al., 2020; Soto & Schurr, 2021). These models convert seismic waveforms into time series known as image functions (or characteristic functions), which quantify the probability of a P-phase, S-phase, or noise being present in the waveform (Zhu & Beroza, 2019). An example image function generated from single-station waveforms is shown in Figure 6 in Section 5.1. A phase detection is declared if the probability is above a certain threshold.

We select the Generalised Phase Detection (GPD) model, which was trained to detect P- and S-phases from events with magnitudes between −0.81 and 5.7 and epicentral distances of less than 100 km (Ross et al., 2018). The model was trained using hand-labelled data archives from the Southern California Seismic Network (California Institute of Technology and United States Geological Survey Pasadena, 1926), which included 1.5 million P- and S-wave seismograms and an equal number of 4-second noise windows.

Step 2: Stacking and migration – Qseek

To associate the neural network’s phase annotations with a coherent seismic source, we apply the stacking and migration approach implemented in the Qseek software (Isken et al., 2025).

Figure 3 illustrates the stacking and migration approach, where the seismograms are annotated with the first arrivals of the P- and S-phases using a pre-trained neural network. The magnitude of each annotation represents the model’s certainty in the phase arrival. These annotations are back-projected onto a subsurface grid according to theoretical travel times and stacked at each node. The node with the most constructive stack corresponds to the maximum semblance (i.e. coherence of arrivals). If this value exceeds a predefined threshold, the algorithm indicates a detection and assigns the corresponding node as the most likely location of the seismic event.

Figure 3. Conceptual 2D illustration of the stacking-and-migration procedure used in Qseek (Isken et al., 2025). (A.1) Image functions recorded at the receiver array (yellow triangles) after shifting each trace by the predicted P-wave travel time from the blue candidate node shown in B. The pulses remain misaligned, and the resulting stack (C.1) displays low coherence. (A.2) Image functions shifted by the travel times from the red candidate node in B. The shifts align the pulses coherently, producing a sharp, high-amplitude peak in the corresponding stack (C.2). B. Two-dimensional search grid beneath the array. The colour map denotes the maximum normalised semblance obtained at each node. The blue and red squares indicate the two candidate nodes evaluated in panels (A.1) and (A.2).

We describe the selection of preprocessing parameters and detection thresholds in Section 5.1.

3.3 Hypocentre estimation

Hypocentre inversion aims to estimate the spatial location x_s (often-times together with the origin time T⁰) at which energy is released during a seismic event using the arrival times of seismic phases at various seismic stations. It involves solving both the forward and the inverse problem. In the forward problem, theoretical seismic-phase arrival times are computed, simulating the travel times of seismic waves originating from a (potential) hypocentre to various receiver stations. The inverse problem seeks to identify the solution(s) that minimise the differences between theoretical phase arrivals (t_calc) and observed phase arrivals (t_obs). Note that different hypocentres may satisfy the fitting criterion, and small perturbations in the input data (e.g. arrival times) can lead to significant changes in the inferred hypocentre. Therefore, we are interested not only in the model parameters that best fit the observed data but also in the overall volume of solutions and their associated uncertainties.

Observed phase arrivals and associated uncertainty

An observed phase arrival corresponds to the (picked) onset time of the first measurable seismic energy recorded on a seismogram. Due to potential errors in picking the onset time, it is more accurate to represent a phase arrival as a probability density function (PDF) rather than a scalar value. Following Lomax et al. (2009), we use a normal distribution to describe an observed phase arrival. The key parameters for the hypocentre estimation are the mean of the distribution, t_obs, and the corresponding standard deviation, σ_obs. Figure 4 illustrates how these terms are derived from the picking process and the corresponding observed phase-arrival distribution. We manually reviewed and selected the phase arrival times and their associated uncertainty.

Figure 4. Illustration of a p-phase arrival time and associated picking uncertainty. The blue line represents a seismic trace, while the green line is a normal distribution that describes the uncertainty of the phase arrival time. Note that the normal distribution is scaled for visual clarity, and its integral equals one. The label on the vertical axis only applies to the seismic trace amplitude.

Theoretical phase arrival times

Estimating the hypocentre requires a model to compute theoretical phase arrival times, which are then compared with the observed phase arrival times. The theoretical phase arrival, t_calc, results from the solution of the forward problem. To compute t_calc, we require the receiver’s location, a hypocentre location (our guess), and the forward function u describing the propagation time from the hypocentre’s location to the receiver’s location. The theoretical phase arrival time is defined as:

In this study, we adopt an infinite frequency approximation (Lin & Ritzwoller, 2011). Specifically, u exploits the fast Marching Method (FMM) (Sethian & Popovici, 1999), implemented in the PyKonal Python Package (White et al., 2020). Besides being contingent on the validity of the infinite frequency approximation, the accuracy of the estimated travel time depends strongly on the accuracy of the velocity model v_model.

Travel-time uncertainty

The accuracy of the velocity model v_model directly affects the accuracy of the theoretical phase arrival times, which in turn impacts the reliability of the hypocentre inversion. The velocity model should closely represent accurate velocities in the subsurface and, as such, correctly account for the interaction of the phases with the different lithostratigraphic units and their geometry. Variations in subsurface properties, such as layer thicknesses or velocity gradients, can lead to deviations between the theoretical and observed phase arrivals. To account for these variations, the velocity model must incorporate both mean velocity values and measures of variability, such as standard deviations or probability distributions. These can be used as inputs to estimate the uncertainties of the theoretical phase arrivals.

The following steps detail the procedure we use to estimate the standard deviation of σ_calc:

1. Retrieve the velocity profile: Select a 1D velocity profile with well-defined mean velocity values and standard deviations for each geological unit (see Section 4.2). For this study, the profile is taken from the surface location of Kwintsheul’s geothermal doublet.

2. Monte Carlo Sampling: Use the mean velocity values and their standard deviations to generate multiple realisations of the velocity model, perturbing the mean velocities by sampling from their respective distributions.

3. Recalculate theoretical phase arrivals: For each perturbed velocity model, recompute the t^calc at the maximum likelihood hypocentre.

4. Quantify variability: Compute the standard deviation of the recalculated t_calc values across all realisations. This standard deviation represents σ_calc.

The resulting σ_calc quantifies the uncertainty in the theoretical phase arrival times caused by variability in the velocity model. This value is subsequently used in the likelihood function to account for travel-time uncertainties in hypocentre estimation. Note that we did not consider potential correlations of the calculated travel times between different stations. Ideally, this would have been the case, but to limit the computational burden, we chose not to do so.

Equal-differential-time likelihood function

A likelihood function is used to quantify how likely it is that a specific hypocentre could explain the observed phase arrival times. The Equal-Differential-Time formulation (Font et al., 2004) is an approximation that allows removing the origin time from the unknown parameters. It is also robust in case of outliers (Lomax et al., 2009). If a hypocentre x_s is perfectly determined within the Earth model, and in the absence of noise, the time difference between the calculated phase arrivals at two seismic stations a and b should be equal to the difference between their corresponding observed phase arrivals, i.e.:

A reliable hypocentre estimation requires accounting for the uncertainties related to the theoretical phase arrivals and observed phase arrivals. Accounting for these terms, the likelihood of the hypocentre becomes:

where L(x_s) estimates the Gaussian likelihood of the hypocentre x_s. The terms and are the observed arrival times at seismic stations a and b, while and are their corresponding modelled arrival times. The variances of the observed and theoretical phase arrivals of stations a and b are denoted by and , respectively.

For stations a and b, the total uncertainty is expressed as:

where and represent the uncertainties in the observed phase arrivals at stations a and b, respectively, and and correspond to the uncertainties in the theoretical phase arrivals. These combined uncertainties are used in calculating the likelihood function L(x_s), which measures the agreement between the observed and theoretical arrival times for a given hypocentre x_s.

L(x_s) reaches its maximum value when the differentials of the observed and theoretical phase arrivals are equal, hence the term Equal Differential Time (EDT). Since the summation over observations occurs outside the exponential, the EDT PDF attains its highest values at locations where the most observation pairs are satisfied, making it robust against outliers (Lomax et al., 2009). In addition, the EDT PDF is independent of the earthquake’s origin time.

Complete probabilistic solution – posterior sampling

Tarantola and Valette (1982) introduced a general methodology to obtain a complete probabilistic solution of a seismic event’s hypocentre by computing the posterior probability distribution (PPD):

Here, m is the model parameter vector, which includes the hypocentre’s spatial coordinates (x, y, z). In this equation, κ is a normalisation constant that ensures the posterior integrates to one, ρ(m) is the prior PDF representing prior knowledge of the model parameters, and L(m) is the likelihood function representing how well the observed data fits a given set of model parameters. In our case, the likelihood is estimated through equation (3), which effectively removes the origin time from m (as explained in Section 3.3). It should be understood that this likelihood is an approximation of the true Gaussian likelihood function in the sense that the correlation between different arrival time differences is ignored. In principle, these can be accounted for using a covariance matrix (Spetzler et al., 2024). This, however, would require us to limit the number of arrival-time differences to N − 1 (where N is the number of stations for which a P-wave arrival-time pick is available). In this study, we choose to adopt the likelihood proposed by Lomax et al. (2009), which is based on all arrival-time differences (N (N − 1)/2). In principle, any set of model parameters can be arbitrarily chosen to calculate theoretical phase arrivals. However, prior information, often based on fundamental laws or physical constraints, can indicate whether a set of parameters is feasible. Here, we use a uniform prior probability density that assumes that any model parameter has an equal probability of explaining the observed data:

That is, we assign equal a priori probabilities to equal volumes. Note that the prior is not completely uninformative, as we still choose a search area.

The spatial domain must be discretised into a grid to perform numerical computations. Each grid cell represents a surrounding region (Δx Δy Δz), which must be sufficiently small to ensure the probability distribution remains approximately constant within the cell (Mosegaard & Tarantola, 1995). Consequently, the probability associated with a point in the grid denotes the probability that the hypocentre falls within the surrounding region represented by the grid cell. The marginal PDF of the hypocentre can be estimated by integrating over specific parameters in the model parameter space. The marginal PDF of the epicentral location can be obtained by integrating over the depth parameter z as:

where the limits of integration represent the plausible range of hypocentral depths.

The marginal PDF for the depth, π_z(z), can be obtained by integrating over the epicentral coordinates x and y from the full posterior distribution. This process effectively reduces the three-dimensional posterior PDF to a one-dimensional distribution along the depth axis, which allows us to understand the uncertainty specifically related to the depth of the hypocentre. The depth PDF is computed as:

where the integration bounds represent the plausible spatial extent of the epicentre.

The spatial uncertainty of the seismic events is quantified by analysing the depth’s marginal PDF π_z. We define the uncertainty range as the interval where the likelihood values exceed 95% of the maximum likelihood, ensuring that it encompasses the most probable depth values. We refer to this range as the 95% confidence interval.

To estimate this range, we follow these steps:

1. Identify the maximum likelihood value, L^max, from the 3D posterior distribution.

2. Define a threshold at 95% of L^max, above which depth solutions are considered significant.

3. Calculate the mean likelihood over x and y. The depth uncertainty range is determined as the interval [z_min, z_max] where the mean likelihood exceeds the threshold.

4. The depth uncertainty is given by Δz = z^max − z^min.

The probabilistic hypocentre determination provides the spatial coordinates and depth of each event, along with their associated uncertainties. These results serve as essential inputs for the next step: estimating event magnitudes.

3.4 Local-magnitude estimation

We estimate the local magnitudes (M_L) following the methodology described by Dost et al. (2004). This approach uses the observed peak amplitudes from seismic waveforms, corrected for distance-dependent attenuation. The attenuation correction factor was calibrated using a set of induced earthquakes at approximately 3 km depth in the Groningen area (Dost et al., 2004).

The procedure for estimating M_L involves the following steps:

1. Waveform preparation: Seismic waveforms are detrended and preprocessed using instrument-response corrections to obtain Wood-Anderson-equivalent amplitudes. We use only horizontal components of the waveform, corresponding to S-wave arrivals, for magnitude estimation.

2. Peak amplitude extraction: For each horizontal trace, the absolute peak amplitude (A_WA) is identified. Here, A_WA is the maximum averaged horizontal-displacement amplitude of a simulated Wood-Anderson instrument, expressed in millimetres.

3. Local-magnitude calculation: The local magnitude for each station is computed as:

   

   (Dost et al., 2004), where A_WA is the observed peak amplitude of the simulated Wood-Anderson displacement, and R is the hypocentral-distance correction term.

4. Station filtering: Magnitudes are calculated for all stations that recorded the event, and the median magnitude is used as the event’s M_L. Stations with anomalous deviations from the event-wise mean are excluded.

4. Data

4.1 Waveform data

The temporary seismic network in Kwintsheul was operational from 22 July 2019 to 9 November 2019 (Muntendam-Bos et al., 2022; Naranjo et al., 2022). The network consisted of 30 three-component force-balance Seismotech geophone sensors, which recorded at a bandwidth of 0.2–100 Hz at 250 sps. These were installed on the surface. The geometry of the network consisted of two intersecting lines, each composed of 13 stations, covering an area of approximately 3.8 km². These stations were installed with an average in-line spacing of 150 meters. The layout was designed to record ambient seismic noise for consecutive illumination analysis and application of body-wave seismic interferometry (e.g. Panea et al., 2014). Additionally, an outer ring of four peripheral stations surrounded the array, encompassing an 18 km² area around the injection point of the geothermal doublet. The goal of the peripheral stations was to increase the azimuthal coverage for the location of events and a more robust estimation of hypocentral depths. The network’s layout and the location of the geothermal doublet are shown in the inset of Figure 1.

As explained in Section 3.1, we assess the completeness of the available seismic waveform data. Figure 5 provides an overview of the data, with vertical lines indicating gaps in the data that directly affect the detection results and hypocentre estimations (see Section 5).

Figure 5. An overview of seismic waveform data visualised using Snuffler (Heimann et al., 2017). Continuous gray colours indicate sections without gaps, while vertical lines indicate interruptions in the data.

The most significant data gaps are observed at stations 027, 029, and 030. These stations are part of the outer ring designed to improve azimuthal coverage and extend the array’s reach (Figure 1). Specific details regarding these gaps include:

• Station 029 stopped recording on 18 July 2019.
• Station 027 stopped recording on 16 August 2019.
• Station 030 experienced intermittent gaps and stopped recording on 29 September 2019.

The effects of these gaps on detection capabilities and hypocentre estimations are further detailed in Section 5.

4.2 Seismic-velocity model

We compute theoretical phase arrivals and travel-time uncertainties (Sections 3.3 and 3.3) using Velmod 3.1 (Pluymaekers et al., 2017). Velmod provides both the mean velocities and corresponding uncertainties for each geological unit. It is parameterised as follows:

where is the P-wave velocity at depth z, is the mean velocity at the top of the unit, is the standard deviation of the velocity, and k is the velocity gradient within the unit. For additional details on retrieving the velocity model, refer to Section 3.1.

5. Results

5.1 Event detections

The magnitude Md 0.16 event reported in Muntendam-Bos et al. (2022) raised the question of whether more seismicity occurred below Kwintsheul. We, therefore, applied the TM detection routine (see Section 3.2) using the waveforms from the reported event as templates. To enhance the signal, we apply a band-pass filter between 4 Hz and 20 Hz and use templates with a length of 1 s following the P- and S-wave phase arrival time. We use templates for each station channel (i.e. P-wave for vertical components and S-wave for horizontal components). As a result, we identified five additional seismic events with near-identical waveforms. Together with the event reported in Muntendam-Bos et al. (2022), this gives six events.

Building on the six identified seismic events, we calibrate the machine-learning detection routine. Among several deep-learning arrival-time picking models, we select the Generalized Phase Detection (GPD) model (Ross et al., 2018), whose training dataset includes low-magnitude events (see Section 3.2 for details). To ensure consistency with the features learned by the GPD model, we apply the same preprocessing: a high-pass filter at 2 Hz and resampling to 100 Hz. Note that these processing parameters deviate from those of the TM routine. We select P- and S-pick thresholds of 0.3 based on the analysis of image functions for the identified events. Figure 6 shows an example image function generated using these threshold values. When the image function exceeds the threshold, a phase arrival is annotated in the time series.

Figure 6. Seismic waveforms of the (A) East, (B) North, and (C) vertical components from station 001 of the temporary seismic array installed in Kwintsheul. (D, E) Corresponding annotated image functions of P- and S-waves using the GPD model with P- and S-pick thresholds of 0.3.

The annotated image functions are the input for the stacking and migration approach explained in Section 3.2. To improve the accuracy of stacking and migration, we define a 1D velocity model retrieved at the location of the injection well. This velocity model is used by Qseek to estimate the travel-time shifts (see Figure 3A). We select a maximum semblance threshold of 0.8, which enables us to detect five of the six events used for calibration.

In total, the machine-learning detection routine identified 65 events that exhibit coherent seismic signals across the seismic array. We classified these events based on their waveform characteristics and spatial distribution. The detection statistics are summarised in Figure 7, and further details on each category can be found in Section 5.2.

Figure 7. Summary of machine-learning detection results. The left chart categorises the trigger detections based on station coverage. In blue, detections during periods of poor station coverage (<50%); in green, detections during adequate station coverage (>50%). The middle chart highlights the filtered dataset (after discarding the events on poor-station-coverage dates), displaying coherent events (green) versus false (incoherent) detections (blue). The right chart classifies these coherent events by event characteristics (green = microseismic cluster at 2.3 km depth, blue = events to the east, grey = controlled explosion, and red = other events). These groups are explained in Section 5.2.

A limitation of the machine-learning detection routine is the disproportionately large number of detections produced when fewer sensors are recording. After the association step outlined in Section 3.2, the detection workflow yielded 17,108 detections. However, 16,812 occurred on just six specific dates when station coverage was particularly low (<50%). These dates and the number of sensors operating at those times are summarised in Table 1.

Table 1. Summary of days with fewer than 16 operational stations and elevated detection numbers

Date	Number of available stations	False detections
29 June 2019	6	2786
9 July 2019	15	36
18 July 2019	13	119
19 July 2019	4	4386
20 July 2019	5	8409
8 November 2019	7	1076

5.2 Waveform characterisation

The 65 events detected in Kwintsheul can be categorised based on their waveform attributes into (1) microseismic cluster, (2) eastern events, (3) other events, and (4) a controlled explosion.

1. Microseismicity Cluster: Corresponds to six events that exhibit clear, impulsive P- and S-wave arrivals with frequencies between 5 and 50 Hz (see Figure 8) and P- to S-wave delays of approximately 1.8 seconds. The character of the phase arrivals is impulsive, with a clear, sudden onset. P-waves are most prominent on the vertical component recordings, whereas the S-waves are most pronounced on the horizontal components.

2. Eastern Events: Correspond to 55 coherent events coming from the East of the array. These events exhibit frequencies from 1 to 25 Hz for P-waves and 1 to 12 Hz for the later arriving waves (delayed by approximately 4 seconds with respect to the impulsive P waves; see Figure 9A and 9C). In an attempt to estimate their hypocentres, we found that they appear to originate at or close to the Earth’s surface. Since there is no mechanical argument for having earthquakes originate at depths less than approximately 200 meters (unconsolidated sediments do not allow for seismic stress release), we believe that they, in fact, originated at the Earth’s surface. The lack of a clear S-wave arrival supports this conclusion. We did not include the hypocentre inversion results because the high azimuthal gap (>250) renders the posterior very broad (although it is clear that they come from the East).

3. Other Events: Correspond to three events that exhibit three pulses of low-frequency waves in the vertical component. The spectrogram analysis (Figure 9B and 9D) reveals that the energy is predominantly concentrated along the vertical component, with frequencies ranging between 5 and 10 Hz. The energy travels primarily vertically as the signals arrive almost simultaneously at all stations. There is, furthermore, no clear evidence of an S-wave associated with these events.

4. Controlled Explosion: Corresponds to a controlled explosion in the North Sea on 2019-10-15 at 17:47:11, also reported in the International Seismological Centre’s On-Line Bulletin (ISC, 2025) with Event ID 618929881.

Figure 8. Analysis of the seismic event on 2019-07-15T12:11:51. (A) Waveforms of the event bandpass filtered between 4 and 20 Hz. (B, D) Map views and cross-sections showing the probability distribution of the event location are displayed on the Viridis colour scale, with green indicating the maximum likelihood location. (C) Spectrograms representative of other stations, exemplified here for station 023, show the vertical, north, and east components. The colour scale represents the logarithm of the power spectral density, with red indicating higher values.

 

Figure 9. Analysis of eastern events and non-seismic events. (A) Waveforms of one of the eastern events recorded on 2019-08-03T04:54:18. (B) Spectrogram of station 011 showing the vertical, north, and east components of the eastern event. (C) Waveforms of one of the non-seismic events detected recorded on 2019-10-16 at 21:26:40. (D). Spectrogram of station 023 showing the vertical, north, and east components of the non-seismic event.

5.3 Hypocentres and magnitudes

To investigate the spatial distribution of the microseismicity cluster, we derived their maximum-likelihood hypocentre solutions with 95% confidence intervals as error bars shown in Figure 10A–C. The hypocentres of the microseismicity cluster are located at approximately 2.3 km depth, with vertical uncertainties ranging from 267 to 735 m (Table 2).

Figure 10. Hypocentres that belong to the microseismic cluster. (A, B, C) Map showing the spatial distribution of event hypocentres with associated 95% confidence interval shown as error bars. The background image corresponds to the 3D velocity model used for the inversion, and its corresponding colour scale is shown in the legend on the right side. (D) Projection in 2D showing the localized events, previously mapped faults (S. Peeters, personal communication, November 2024), and the seismic stations. The colour scale indicates the depth at which the events and faults are located.

Table 2. Characterization of the detected microseismic events. Origin times are given in Coordinated Universal Time (UTC). Hypocentre coordinates are reported in the Dutch national reference system (EPSG: 28992). M_Ldenotes the local magnitude, and Var(M_L) represents the variance of the magnitude estimates across multiple stations. The parameters σ_x, σ_y, and σ_z correspond to the 95% confidence intervals of the hypocentre location in each coordinate direction.

Date	Time	x [m]	y [m]	z [m]	Az. gap	ML	Var(ML)	σx [m]	σy [m]	σz [m]
2019-06-23	08:02:08	79483.25	447004.24	–2314.38	93	–2.04	0.13	602	363	267
2019-07-14	08:48:31	79483.30	446967.93	–2354.52	143	–1.60	0.21	662	399	294
2019-07-15	12:11:51	79543.51	446986.10	–2314.38	95	–1.79	0.18	632	381	280
2019-09-11	08:51:06	79212.04	446986.08	–2274.25	149	–2.00	0.13	1627	980	722
2019-10-02	11:46:50	79392.84	447004.24	–2100.33	150	–2.11	0.13	1446	871	642
2019-10-06	01:20:01	79392.90	446986.08	–2153.85	170	–2.52	0.18	1657	998	735

Events on 2019-09-11, 2019-10-02 and 2019-10-06 have different depth hypocentre values (~2.1 km), as shown in Table 2. These events were not recorded by any peripheral station, which increased their azimuthal gap and, consequently, their depth uncertainty, as shown in Appendix A.

To contextualise the hypocentre solutions with the known faults in the Kwintsheul area, we compare them in Figure 10D against the mapped faults, the geothermal doublet, and the seismic stations. The interpretation of these mapped faults is derived from a regional 3D seismic dataset reprocessed in 2012, covering approximately 1200 km² in the WNB (Merrifield, 2012), and was provided by S. Peeters (personal communication, November 2024). The six seismic events lie close to a local fault as well as the injection well, as shown in Figure 10D. We estimate the minimum distance from each event to the nearest point on the closest fault, finding values ranging from 265 to 434 m. However, since the uncertainties in the X, Y, and Z coordinates are within this range, it is not possible to determine whether the events originated on the fault.

The local magnitudes of the microseismicity cluster range from –1.60 to –2.52, computed with the site-specific local-magnitude formula given in Equation 9. The deviation of station-specific M_L’s with respect to the median varies from 0.13 to 0.18 M_L, which shows a relatively low variability among the different stations. This is not surprising, given that the stations are relatively close to each other.

6. Discussion

A primary motivation for this study was to develop a workflow to monitor low-magnitude seismic events in the Netherlands, particularly those around areas with geothermal operations. Although the nationwide KNMI network effectively detects events above the regional and spatially varying magnitude of completeness (Ruigrok et al., 2023), lower magnitude events will go unnoticed. Detecting and accurately locating these low-magnitude events is crucial for understanding seismicity patterns, informing network designs, and ensuring the safe expansion of geothermal operations. Our results confirm that such low-magnitude seismic events occur in the study area.

We demonstrated that local dense seismic arrays are essential for microseismic-event detectability. False detections increased sharply when the array was partially or poorly configured, highlighting the importance of robust station coverage. Moreover, noisy urban environments benefit from a multi-method approach that includes deep-learning pickers and TM (Panebianco et al., 2023; Sugan et al., 2023; Diaferia et al., 2024). Using the multimethod approach introduced in Section 3.2, we identified 65 events that exhibited coherent signals during the 4-month monitoring period.

We found that the seismic event, reported in Muntendam-Bos et al. (2022), was not an isolated occurrence but part of a repetitive microseismic sequence. This follows from applying the TM method, which identified five similar events. The nearly identical waveform patterns indicate a common source region and mechanism. The six seismic events occur at a depth of approximately 2.3 km near both a mapped fault system and the injection well KW-GT-01. The depth uncertainties, however, prevent a definitive interpretation of whether these events are fault-related or induced by the injection. The limited duration of the monitoring period further complicates this distinction. Longer-term monitoring would help clarify whether these events are truly linked to seismic stress release in the geothermal reservoir, activities at the well, or natural tectonic activity.

In addition, we detected seismicity originating to the east of the array, referred to as Eastern Events (or Cluster 2 events). These events exhibit impulsive P-wave arrivals, but they lack an impulsive S-wave. The attenuation of high-frequency S-waves suggests a more distant source for Cluster 2 events. The high azimuthal gap in the current array configuration prevented us from properly estimating their hypocentres, but an origin at the Earth’s surface is likely (though not confirmed). Clearly, low detection thresholds are needed to identify low-magnitude seismicity. These low thresholds result in a higher number of false detections, which requires manual inspection. In regions with high seismicity rates, such manual inspection may be infeasible. Therefore, we recommend retraining the machine-learning pickers on site-specific datasets to improve pick accuracy, enable lower thresholds, and, thereby, reduce false detections. This strategy is particularly suitable where permanent seismic stations and high seismicity rates are available.

Our findings highlight a critical limitation of the current array design: excessive reliance on the outer ring of stations. Although the two lines of aligned sensors provide redundancy in signal detection, they do little to improve azimuthal coverage or depth resolution. This array design is, therefore, suboptimal for passive seismic studies that require accurate event locations (at reservoir depths of two or three kilometres). An optimized network targeting the already localized seismic events is recommended for future campaigns (e.g. Maurer et al., 2010; Toledo et al., 2020; Esquivel-Mendiola et al., 2022). For large-scale operations, including borehole geophones would aid in estimating the depths of hypocentres.

Beyond the need for better station coverage, our study underscores the importance of refining the velocity model used for event location. Although the publicly available P-wave velocity model (Velmod 3.1) was instrumental for this study, an S-wave velocity model is needed to enhance the depth resolution and reduce location uncertainties (e.g. Spetzler et al., 2024).

Our analysis also highlights the challenges of magnitude estimation for low-magnitude events. The Groningen-calibrated magnitude formula (Equation 9), used here as a first-order approximation, may result in systematic under- or over-estimation of event magnitudes in South Holland, where the subsurface properties are different. In addition, Equation 9 is calibrated using recordings from borehole geophones at 200 meters depth. It should, in principle, be used for particle motions recorded at that (or close to that) depth. Establishing a specialised local magnitude scale for this region would require a dedicated network design, knowledge of site-specific attenuation properties, and additional data collected over a longer monitoring period (i.e. more earthquakes).

Taken together, our findings address an important gap between large-scale seismic detection frameworks and the finer resolution needed for local hazard assessment. Moving forward, efforts to refine velocity models, improve monitoring networks, and establish regional magnitude scales will enhance our ability to capture and interpret low-magnitude seismicity.

7. Conclusions

We presented a comprehensive workflow for monitoring seismicity, specifically designed to detect and characterise low-magnitude events in urban areas. We addressed key challenges in seismic monitoring by integrating data preparation, template-matching detection, machine-learning-based detection, and probabilistic hypocentre estimation. Applying this workflow to the Kwintsheul area, we detected six seismic events near a local fault and near the bottom of injection well KW-GT-01, although spatial uncertainties remain in the order of hundreds of meters. In addition, we detected 59 events with coherent signals, but inadequate azimuthal coverage hindered their accurate characterisation. Our findings highlight the need for more spatially distributed networks, refined local velocity models, and region-specific magnitude calibrations to enhance the accuracy and reliability of seismic monitoring. Our workflow provides a scalable and adaptable solution for improving seismic monitoring in urban environments.

Acknowledgements

We want to thank Paul van Schie, Mathieu de Bas, and the Nature's Heat project for encouraging us to monitor the project and helping us install and maintain the network. We thank Stefan Peeters (TNO Geological Survey of the Netherlands) for providing the fault interpretation used for Figure 10. DN has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska Curie grant agreement No 956965. Map data is copyrighted by OpenStreetMap contributors (2017). Finally, we would like to express our sincere gratitude to the editor and anonymous reviewers for their constructive comments and suggestions, which significantly improved the quality of the manuscript.

References

Boersma, Q.D., Bruna, P.O., de Hoop, S., Vinci, F., Moradi Tehrani, A. & Bertotti, G., 2021. The impact of natural fractures on heat extraction from tight Triassic sandstones in the West Netherlands Basin: a case study combining well, seismic and numerical data. Netherlands Journal of Geosciences 100: e6. DOI: 10.1017/njg.2020.21.

California Institute of Technology and United States Geological Survey Pasadena, 1926. Southern California seismic network. https://www.fdsn.org/networks/detail/CI.

Chamberlain, C.J., Hopp, C.J., Boese, C.M., Warren-Smith, E., Chambers, D., Chu, S.X., Michailos, K. & Townend, J., 2017. EQcorrscan: repeating and near-repeating earthquake detection and analysis in python. Seismological Research Letters 89(1): 173–181. DOI: 10.1785/0220170151.

Diaferia, G., Valoroso, L., Improta, L. & Piccinini, D., 2024. A high-resolution seismic catalog for the southern Apennines (Italy) built through template-matching. Geochemistry, Geophysics, Geosystems 25(3): e2023GC011160. DOI: 10.1029/2023GC011160.

Dost, B., van Eck, T. & Haak, H.W., 2004. Scaling of peak ground acceleration and peak ground velocity recorded in the Netherlands. Bollettino Di Geofisica Teorica Ed Applicata 45: 153–168.

Duin, E., Doornenbal, J., Rijkers, R., Verbeek, J. & Wong, T., 2006. Subsurface structure of the Netherlands – results of recent onshore and offshore mapping. Netherlands Journal of Geosciences – Geologie en Mijnbouw 85(4): 245–276. DOI: 10.1017/S0016774600023064.

Ellsworth, W., 1995. Characteristic earthquakes and long-term earthquake forecasts: implications of central California seismicity. In: Cheng, F. & Sheu, M.S. (eds.): Urban disaster mitigation: the role of engineering and technology. Pergamon (Oxford): 1–14. DOI: 10.1016/B978-008041920-6/50007-5

Esquivel-Mendiola, L.I., Calò, M., Tramelli, A. & Figueroa-Soto, A., 2022. Optimization of local scale seismic networks applied to geothermal fields. The case of the acoculco caldera, Mexico. Journal of South American Earth Sciences 119: 103995. DOI: 10.1016/j.jsames.2022.103995.

Font, Y., Kao, H., Lallemand, S., Liu, C.-S. & Chiao, L.-Y., 2004. Hypocentre determination offshore of eastern Taiwan using the maximum intersection method. Geophysical Journal International 158(2): 655–675. DOI: 10.1111/j.1365-246X.2004.02317.x.

Groos, J.C. & Ritter, J.R.R., 2009. Time domain classification and quantification of seismic noise in an urban environment. Geophysical Journal International 179(2): 1213–1231. DOI: 10.1111/j.1365-246X.2009.04343.x.

Heimann, S., Kriegerowski, M., Isken, M., Cesca, S., Daout, S., Grigoli, F., et al. (2017). Pyrocko – an open-source seismology toolbox and library [Software]. GFZ Data Services, Potsdam. https://doi.org/10.5880/GFZ.2.1.2017.001.

Heimann, S., Kriegerowski, M., Isken, M., Cesca, S., Daout, S., Grigoli, F., Juretzek, C., Megies, T., Nooshiri, N., Steinberg, A., Sudhaus, H., Vasyura-Bathke, H., Willey, T. & Dahm, T., 2017. Pyrocko – an open-source seismology toolbox and library.

Houtgast, R. & van Balen, R., 2000. Neotectonics of the Roer Valley Rift System, the Netherlands. Global and Planetary Change 27(1): 131–146. DOI: 10.1016/S0921-8181(01)00063-7.

Hoyer, S. & Hamman, J., 2017. xarray: N-D labeled arrays and datasets in Python. Journal of Open Research Software 5: 10. DOI: 10.5334/jors.148.

ISC, 2025. International Seismological Centre. On-line Bulletin. DOI: 10.31905/D808B830.

Isken, M., Heimann, S., Niemz, P., Münchmeyer, J., Cesca, S., Vasyura-Bathke, H. & Dahm, T., 2025. Qseek: a data-driven framework for automated earthquake detection, localization and characterization. Seismica 4(1): 1283. DOI: 10.26443/seismica.v4i1.1283.

KNMI, 1993. Netherlands Seismic and Acoustic Network [Data set]. Royal Netherlands Meteorological Institute (KNMI), De Bilt. DOI: 10.21944/e970fd34-23b9-3411-b366-e4f72877d2c5.

KNMI, 2023. Earthquakes – complete catalogue for the Netherlands and near surrounding [Data set]. Royal Netherlands Meteorological Institute (KNMI), De Bilt. Available at: https://dataplatform.knmi.nl/dataset/aardbevingen-catalogus-1.

Kraaijpoel, D. (2025). TNO/PRESEIS-dgm-velmod-sampler (v2025.1) [Software]. Zenodo. https://doi.org/10.5281/zenodo.14698789.

Lin, F.-C. & Ritzwoller, M.H., 2011. Helmholtz surface wave tomography for isotropic and azimuthally anisotropic structure. Geophysical Journal International 186(3): 1104–1120. DOI: 10.1111/j.1365-246X.2011.05070.x.

Lomax, A., Michelini, A. & Curtis, A., 2009. Earthquake location, direct, global-search methods. Springer, New York: 1–33. DOI: 10.1007/978-3-642-27737-5_150-2.

Maurer, H., Curtis, A. & Boerner, D.E., 2010. Recent advances in optimized geophysical survey design. Geophysics 75(5): 75A177–75A194. DOI: 10.1190/1.3484194.

Menke, W., 1999. Using waveform similarity to constrain earthquake locations. Bulletin of the Seismological Society of America 89(4): 1143–1146. DOI: 10.1785/BSSA0890041143.

Merrifield, T. (2012). R-2823 R-2823 Rotterdam Donkersloot Pre-SDM. Shell UIE, March 2012. Internal confidential Shell UIE report.

Mosegaard, K. & Tarantola, A., 1995. Monte Carlo sampling of solutions to inverse problems. Journal of Geophysical Research: Solid Earth 100(B7): 12431–12447. DOI: 10.1029/94JB03097.

Mousavi, S.M., Ellsworth, W.L., Zhu, W., Chuang, L.Y. & Beroza, G.C., 2020. Earthquake transformer – an attentive deep-learning model for simultaneous earthquake detection and phase picking. Nature Communications 11: 3952. DOI: 10.1038/s41467-020-17591-w.

Muntendam-Bos, A.G., Hoedeman, G., Polychronopoulou, K., Draganov, D., Weemstra, C., Van Der Zee, W., Bakker, R.R. & Roest, H., 2022. An overview of induced seismicity in the Netherlands. Netherlands Journal of Geosciences 101: e1. DOI: 10.1017/njg.2021.14.

Naranjo, D., Draganov, D., Polychronopoulou, K., De Bas, M., & Weemstra, C. (2022). Seismic monitoring of Nature's Heat Geothermal Reservoir in Kwintsheul (Netherlands). Paper presented at the European Geothermal Congress 2022, Berlin, Germany, 17–21 October 2022. https://www.egec.org/wp-content/uploads/2023/02/Naranjo-Microseismic-monitoring-of-Natures-Heat-Geothermal-Project-432_ExtAbstract.pdf

OpenStreetMap Contributors, 2017. Planet dump. https://www.openstreetmap.org.

Panea, I., Draganov, D., Vidal, C.A. & Mocanu, V., 2014. Retrieval of reflections from ambient noise recorded in the Mizil Area, Romania. Geophysics 79(3): Q31–Q42. DOI: 10.1190/geo2013-0292.1.

Panebianco, S., Serlenga, V., Satriano, C., Cavalcante, F. & Stabile, T.A., 2023. Semi-automated template matching and machine-learning based analysis of the August 2020 Castelsaraceno microearthquake sequence (southern Italy). Geomatics, Natural Hazards and Risk 14(1): 2207715. DOI: 10.1080/19475705.2023.2207715.

Pluymaekers, M., Doornenbal, J. & Middelburg, H., 2017. TNO 2017 R11014 with erratum page 67 Final, Velmod-3.1. TNO Report. https://www.nlog.nl/sites/default/files/2018-11/060.26839%20R11014%20with%20erratum%20page%2067%20Doornenbal-final.sec_.pdf.

Ross, Z.E., Meier, M., Hauksson, E. & Heaton, T.H., 2018. Generalized seismic phase detection with deep learning. Bulletin of the Seismological Society of America 108(5A): 2894–2901. DOI: 10.1785/0120180080.

Ruigrok, E., Kruiver, P. & Dos, B., 2023. Construction of earthquake location uncertainty maps for the Netherlands. KNMI number: TR-405. p. 158.

Ruigrok, E., Kruiver, P. & Dost, B. (2023). Construction of earthquake location uncertainty maps for the Netherlands. KNMI Report TR-405, Royal Netherlands Meteorological Institute (KNMI), De Bilt, p. 158.

Sethian, J.A. & Popovici, A.M., 1999. 3-d traveltime computation using the fast marching method. Geophysics 64(2): 516–523. DOI: 10.1190/1.1444558.

Soto, H. & Schurr, B., 2021. DeepPhasePick: a method for detecting and picking seismic phases from local earthquakes based on highly optimized convolutional and recurrent deep neural networks. Geophysical Journal International 227(2): 1268–1294. DOI: 10.1093/gji/ggab266.

Spetzler, J., Ruigrok, E. & Bouwman, D., 2024. Hypocenter uncertainty analysis of induced and tectonic earthquakes in the Netherlands. Journal of Seismology 28: 555–577. DOI: 10.1007/s10950-024-10205-8.

Sugan, M., Peruzza, L., Romano, M.A., Guidarelli, M., Moratto, L., Sandron, D., Plasencia Linares, M.P. & Romanelli, M., 2023. Machine learning versus manual earthquake location workflow: testing LOC-FLOW on an unusually productive microseismic sequence in northeastern Italy. Geomatics, Natural Hazards and Risk 14(1): 2284120. DOI: 10.1080/19475705.2023.2284120.

Tarantola, A. & Valette, B., 1982. Inverse problems = Quest for information. Journal of Geophysics 50: 159–170.

Toledo, T., Jousset, P., Maurer, H. & Krawczyk, C., 2020. Optimized experimental network design for earthquake location problems: applications to geothermal and volcanic field seismic networks. Journal of Volcanology and Geothermal Research 391: 106433. DOI: 10.1016/j.jvolgeores.2018.08.011.

v. Gessel, S., Hintersberger, E., v. Ede, R., ten Veen, J., Doornenbal, H., Diepolder, G., den Dulk, M., Hamiti, S., Vukzaj, N., Çako, R., Prendi, E., Ceroni, M., Mara, A., Barros, R., Tovar, A., Britze, P., Baudin, T., Stück, H., Jähne-Klingberg, F., Jahnke, C., Höding, T., Malz, A., Kristjánsdóttir, S., Þorbergsson, A., Di Manna, P., D’Ambrogi, C., Congi, M., Lazauskiene˙, J., Andriuškevicˇiene˙, G., Baliukevicˇius, A., Jarosin´ski, M., Gogołek, T., Ste˛pien´, U., Krzemin´ska, E., Salwa, S., Habryn, R., Aleksandrowski, P., Szynkaruk, E., Konieczyn´ska, M., Ressurreição, R., Machado, S., Moniz, C., Sampaio, J., Dias, R., Carvalho, J., Fernandes, J., Ramalho, E., Filipe, A., Celarc, B., Atanackov, J., JamŠek Rupnik, P., Shevchenko, A., Melnyk, I. & Lapshyna, A., 2021. The HIKE European fault database (EFDB) compiled in the framework of the Geo-ERA project HIKE (2018–2021). Deliverable GeoERA project HIKE, European Geological Data Infrastructure (EGDI). https://egdi.geology.cz/record/basic/5edf7bd4-9270-4188-b69d-7ddd0a010833.

White, M.C.A., Fang, H., Nakata, N. & Ben-Zion, Y., 2020. PyKonal: a python package for solving the Eikonal equation in spherical and Cartesian coordinates using the fast marching method. Seismological Research Letters 91(4): 2378–2389. DOI: 10.1785/0220190318.

Worum, G., Michon, L., van Balen, R.T., van Wees, J.-D., Cloetingh, S. & Pagnier, H., 2005. Pre-Neogene controls on present-day fault activity in the West Netherlands Basin and Roer Valley Rift System (southern Netherlands): role of variations in fault orientation in a uniform low-stress regime. Quaternary Science Reviews 24(3): 473–488. DOI: 10.1016/j.quascirev.2004.02.020.

Zhu, W. & Beroza, G.C., 2019. Phasenet: a deep-neural-network-based seismic arrival-time picking method. Geophysical Journal International 216(1): 261–273. DOI: 10.1093/gji/ggy423.

Appendix

A Influence of station 030 on hypocentre uncertainty

The influence of Station 030 on hypocentre location accuracy is illustrated in Figure A1, which compares two inversions: one excluding data from Station 030 (Panels A and C) and one including it (Panels B and D). The likelihood probability density function (PDF) is visualised using a Viridis colour scale, where higher values indicate greater likelihood.

Figure A1 Probability density function (PDF) projections for two seismic events. (A and C) show the map view and cross-section of the event on 14 July 2019, while (B and D) show the event on 6 October 2019. Red and blue lines represent the geothermal wells KW-GT01 and KW-GT02. Dark red triangles indicate seismic stations, and dotted lines mark the cross-section positions. PDF contours (Viridis scale) are overlaid on the velocity model used for the inversion.

In the inversion without Station 030 (Figure A1.A and C), the maximum likelihood location is estimated at x = 79355.37, y = 446982.12, and z = −2084.03 m. The uncertainties for this estimate are: azimuthal gap of 168.68◦, depth uncertainty of 470.59 m, x uncertainty of 1060.01 m, and y uncertainty of 638.47 m.

Including Station 030 in the inversion (Figure A1.B and D) significantly reduces uncertainties and shifts the maximum likelihood location. The new estimate is at x = 79506.80, y = 446936.51, and z = −2386.55 m. The uncertainties for this solution are: azimuthal gap of 146.07◦, depth uncertainty of 201.68 m, x uncertainty of 454.29 m, and y uncertainty of 273.63 m.