The Hamiltonian Monte Carlo algorithm is known to be highly efficient when sampling high-dimensional model spaces due to Hamilton's equations guiding the sampling process. For weakly non-linear problems, linearizing the forward problem enhances this efficiency. This study integrates this linearization with geological prior knowledge for optimal results. We test this approach to estimate the source parameters of a 3.4 magnitude induced event that originated in the Groningen gas field in 2019. The source parameters are the event's centroid (three components), its moment tensor (six components), and its origin time. In terms of prior knowledge, we tested two sets of centroid priors. The first set exploits the known fault geometry of the Groningen gas field, whereas the second set is generated by placing initial centroid priors on a uniform horizontal grid at a depth of 3 km (the approximate depth of the gas reservoir). As for the forward problem linearization, we use an approach in which the linearization is run iteratively in tandem with updates of the centroid prior. We demonstrate that, in the absence of a sufficiently accurate initial centroid prior, the linearization of the forward model necessitates multiple initial centroid priors. Eventually, both prior sets yield similar posteriors. Most importantly, however, they agree with the geological knowledge of the area: the posterior peaks for model vectors containing a centroid near a major fault and a moment tensor that corresponds to normal faulting along a plane with a strike almost aligning with that of the major fault.

@en

Carbon capture and storage (CCS) technology is essential to European decarbonisation efforts, and several offshore CO₂ storage projects are being developed in the North Sea. Understanding the geomechanical response to CO₂ injection is key to both the pre-characterisation and operation of a storage reservoir. A thorough assessment of seismicity gives critical insights into the stress field and faulting around reservoirs, both key controls on the geomechanical response to injection. Seismicity also illuminates potential hydraulic pathways for leakage, be it directly by revealing the extent of faults, or indirectly through fractures imaged by measurements of seismic anisotropy. High quality seismicity data is critical to underpin all of these methods of analysis. This paper presents the most complete catalogue of seismicity in the North Sea to date. The combined data are enabling revised assessments of seismic hazard and leakage risk in the North Sea, as well as a better understanding of faulting and stress. This study shows the value of unifying disparate seismicity data, allowing for more accurate seismological analyses. These lay the foundation for better management of risks for not only geologic CO₂ storage, but other offshore industries and infrastructure.

@en

Ocean-bottom seismometers (OBSs) are equipped with seismic sensors that record acoustic and seismic events at the seafloor, which makes them suitable for investigating tectonic structures capable of generating earthquakes offshore. One critical parameter to obtain accurate earthquake locations is the absolute time of the incoming seismic signals recorded by the OBSs. It is, however, not possible to synchronize the internal clocks of the OBSs with a known reference time, given that GNSS signals are unable to reach the instrument at the sea bottom. To address this issue, here we introduce a new method to synchronize the clocks of large-scale OBS deployments. Our approach relies on the theoretical time-symmetry of time-lapse (averaged) crosscorrelations of ambient seismic noise. Deviations from symmetry are attributed to clock errors. This implies that the recovered clock errors will be obscured by lapse crosscorrelations' deviations from symmetry that are not due to clock errors. Non-uniform surface wave illumination patterns are arguably the most notable source which breaks the time symmetry. Using field data, we demonstrate that the adverse effects of non-uniform illumination patterns on the recovered clock errors can be mitigated by means of a weighted least-squares inversion that is based on station-station distances. In addition, our methodology permits the recovery of timing errors at the time of deployment of the OBSs. This error can be attributed to either: i) a wrong initial time synchronization of the OBS or ii) a timing error induced by changing temperature and pressure conditions while the OBS is sunk to the ocean floor. The methodology is implemented in an open-source Python package named OCloC, and we applied it to the OBS recordings acquired in the context of the IMAGE project in and around Reykjanes, Iceland. As expected, most OBSs suffered from clock drift. Surprisingly, we found incurred timing errors at the time of deployment for most of the OBSs.@en

With the development of several CO₂ storage operations in the North Sea, there is a clear need to better characterise the seismic hazard and stress state in the region. Faults and associated fracture sets can act as hydraulic pathways for unintended CO₂ migration, ill-defined stress states can lead to numerous operational difficulties, and induced seismicity will be a clear risk as CO₂ is injected into subsurface reservoirs. Seismicity can reveal the location and extent of faults and fractures, and can be used to invert for the state of stress. Both operators and regulators therefore need a clear understanding of the rate of natural seismicity, to identify and distinguish induced events from natural, and to assess the likelihood of induced fault reactivation. This requires a dedicated, site specific background monitoring programme, as well as a high-quality seismic catalogue for the region around any CO₂ storage operation. Our study has produced the first dedicated seismic catalogue of the North Sea, based on all available data from each of the relevant seismological agencies. This dataset fosters further studies into seismic hazard, leakage risk, and stress state in a region that will be vital for European CO₂ storage efforts in the coming decades.@en

We use the Hamiltonian Monte Carlo (HMC) algorithm to estimate the posterior probability distribution of a number of earthquake source parameters. This distribution describes the probability of these parameters attaining a specific set of values. The efficiency of the HMC algorithm, however, can be improved through the formulation of a geologically constrained prior probability distribution. The primary objective of the presented study is, therefore, to assess the role of the prior probability in the application of the HMC algorithm to recordings of induced seismic events in the Groningen gas field.

@en

Ambient noise seismic tomography has proven to be an effective tool for subsurface imaging, particularly in volcanic regions such as the Reykjanes Peninsula (RP), SW Iceland, where ambient seismic noise is ideal with isotropic illumination. The primary purpose of this study is to obtain a reliable shear wave velocity model of the RP, to get a better understanding of the subsurface structure of the RP and how it relates to other geoscientific results. This is the first tomographic model of the RP which is based on both on- and off-shore seismic stations. We use the ambient seismic noise data and apply a novel algorithm called one-step 3-D transdimensional tomography. The main geological structures in the study area (i.e. covered by seismic stations) are the four NE-SW trending volcanic systems, orientated highly oblique to the plate spreading on the RP. These are from west to east; Reykjanes, Eldvörp-Svartsengi, Fagradalsfjall and Krýsuvík, of which all except Fagradalsfjall host a known high-temperature geothermal field. Using surface waves retrieved from ambient noise recordings, we recovered a 3-D model of shear wave velocity. We observe low-velocity anomalies below these known high-temperature fields. The observed low-velocity anomalies below Reykjanes and Eldvörp-Svartsengi are significant but relatively small. The low-velocity anomaly observed below Krýsuvík is both larger and stronger, oriented near-perpendicular to the volcanic system, and coinciding well with a previously found low-resistivity anomaly. A low-velocity anomaly in the depth range of 5-8 km extends horizontally along the whole RP, but below the high-temperature fields, the onset of the velocity decrease is shallower, at around 3 km depth. This is in good agreement with the brittle-ductile transition zone on the RP. In considerably greater detail, our results confirm previous tomographic models obtained in the area. This study demonstrates the potential of the entirely data-driven, one-step 3-D transdimensional ambient noise tomography as a routine tomography tool and a complementary seismological tool for geothermal exploration, providing an enhanced understanding of the upper crustal structure of the RP.

@en

The timing of the recordings of ocean-bottom seismometers (OBSs) is critical for accurate earthquake location and Earth model studies. GNSS signals, however, cannot reach OBSs deployed at the ocean bottom. This prevents their clocks from being synchronized with a known reference time. To overcome this, we developed OCloC, a Python package that uses time-lapse cross-correlations of ambient seismic noise to synchronize the recordings of large-scale OBS deployments. By simultaneously quantifying deviations from symmetry of a set of lapse cross-correlations, OCloC recovers the incurred clock errors by means of a least-squares inversion. In fact, because non-uniform noise illumination patterns also break the symmetry of (lapse) cross-correlations, we introduce a distance-based weighted least-squares inversion. This mitigates the adverse effect of the noise illumination on the recovered clock errors. Using noise recordings from the IMAGE project in Reykjanes, Iceland, we demonstrate that OCloC significantly reduces the time and effort needed to detect and correct timing errors in large-scale OBS deployments. In addition, our methodology allows one to evaluate potential timing errors at the time of OBS deployment. These might be caused by incorrect initial synchronization, or by rapidly changing temperature conditions while the OBS is sunk to the sea bottom. Our work advances the use of OBSs for earthquake studies and other applications.@en

We report on a Bayesian (i.e., probabilistic) inversion for the shear-wave velocity structure of the Reykjanes peninsula, SW Iceland. Travel times of Rayleigh waves traversing the peninsula served as input to the probabilistic algorithm. These Rayleigh waves were retrieved through the application of seismic interferometry to yearlong recordings of ambient seismic noise. The Reykjanes peninsula is well placed for this technique because it is surrounded by ocean, which implies a relatively uniform seismic noise illumination; the latter being a condition for accurate interferometric surface wave retrieval. The Bayesian algorithm uses a variable model parametrization by employing Voronoi cells in conjunction with a reversible jump Markov chain Monte Carlo sampler. The algorithm is entirely data-driven, meaning that, contrary to conventional deterministic tomographic inversions, the user does not need to define any regularization or parameterization parameters to solve the inverse problem. The geology in the area of interest is characterized by four NE-SW trending volcanic systems, orientated oblique to the divergent plate boundary cutting across the Reykjanes Peninsula. These are from west to east; Reykjanes, Svartsengi, Fagradalsfjall and Krýsuvík, of which all except Fagradalsfjall host a known high-temperature geothermal field. We observe relatively high shear wave velocity patches close to the Earth’s surface (top two kilometers) at the location of these known high-temperature fields. These high velocity anomalies invert to relatively low shear wave velocities (in comparison to shear wave velocities in the same horizontal plane) at depths greater than 3 km. The latter low-velocity anomalies are relatively small below Reykjanes and Svartsengi. At depths of 5 to 8 km, a low-velocity anomaly extends horizontally below Reykjanes and Svartsengi, correlating relatively well with the inferred brittle-ductile transition below the high-temperature fields at 4-5 km depth. The low-velocity anomaly below Krýsuvík is much larger and coincides with a deep-seated low electrical resistivity anomaly. Horizontally, it coincides with the center of an inflation source at 4–5 km depth. For example, in 2010 this resulted in an uplift exceeding 50 mm/year, but several periods of alternating uplift and subsidence associated with increased seismicity have been observed in Krýsuvík since 2009. Our results both confirm and add details to previous models obtained in the area. Our study demonstrates the potential of Bayesian surface wave inversion as a complementary geophysical tool for geothermal exploration. @en

Coda wave interferometry (CWI) holds promise as a technique for concrete stress monitoring. This is because the coda, which consists of multiply scattered arrivals, is the result of propagation through the medium over large distances. As such, it is sensitive to both minute structural changes and small velocity changes in that medium. Previous studies focusing on concrete have predominantly utilized the time-domain-based stretching technique to measure travel-time changes. There is, however, a lack of consensus on how to quantify these changes effectively. In this study, we conduct a systematic comparison between two techniques, namely the stretching technique and the wavelet cross-spectrum (WCS) technique, for measuring stress-induced velocity changes in a cylindrical concrete sample. Our comparison focuses on two key aspects: (i) stability against cycle skipping and (ii) consistency in retrieving velocity changes. Experimental results reveal that both the WCS technique and the stretching technique yield consistent velocity changes. In terms of stability, it is challenging to determine which technique performs better, due to differences in the mechanisms triggering cycle skipping. However, when considering waves with frequencies ranging from 50 kHz to 80 kHz, both techniques exhibit comparable performance. Based on our findings, we offer the following recommendations for utilizing these CWI techniques in concrete stress monitoring: For the stretching technique, selecting the time window length based on the wave frequency and the expected magnitude of velocity change. For the WCS technique, operating it in the frequency band where spectral decomposition shows sufficiently high energy in the signal and can accommodate the expected magnitude of velocity change.

@en

Geothermal energy is a cleaner and more sustainable source of power, which plays a key role in the transition to a low-carbon economy. Sustainable and safe exploitation of geothermal resources, however, depends on our ability to understand and manage the associated seismic risks. In 2018, Nature's Heat geothermal project began operations in Kwintsheul, Netherlands, aiming to supply heat to 64 hectares of greenhouses. Between July and October 2019, a temporary seismic array was installed to monitor for possible seismic activity at the site. Microseismic moment-tensor inversion is a valuable tool for understanding the mechanics and structure of geothermal reservoirs, and for optimizing their exploitation. It can be challenging, though, to apply this technique when there are high levels of ambient seismic noise, as is often the case in geothermal operations in densely populated areas. In this study, we evaluate the feasibility of inverting the centroid moment tensor of microseismic events in Kwintsheul, using probabilistic moment-tensor inversions. We first test the probabilistic inversion using synthetic recordings of ambient seismic noise, after which we apply the technique to the low-magnitude (Md=0.16) event recorded on July 14, 2019. Our results give insight into the challenges and limitations of applying moment-tensor inversion to low-magnitude events in the context of geothermal operations in the Netherlands.@en

In 2018, a geothermal doublet started operating in Kwinstheul, Netherlands, for supplying heat to 64 hectares of greenhouses corresponding to Nature’s Heat joint initiative. This kind of geothermal operation requires extraction, circulation, and reinjection of fluids at a depth of 2.4 km. The reservoir used for the geothermal operation has shown good hydraulic parameters which allow the circulation of the fluid. Several authors agree that this kind of geothermal operation is unlikely to generate felt seismicity, nevertheless, adequate seismic monitoring is critical to guarantee sustainable and safe use of the subsurface. To monitor the operation of Nature’s Heat project, 30 three-component short-period seismic sensors were installed by Delft University of Technology and Seismotech (Greece). A challenge for seismic monitoring in Kwinstheul is the high levels of seismic noise coming from anthropogenic and operational activities. Despite the high background noise levels, a seismic event of Md 0.16 was recorded on July 14, 2019. To understand the relation of the event and improve the safety of the geothermal operation, we are developing an optimized monitoring scheme. @en

In 2018, the geothermal project Nature's Heat started its operations to supply heat to 64 hectares of greenhouses in Kwintsheul, Netherlands. The operation involves the extraction and reinjection of geothermal fluids at a depth of about 2.4km. Several studies suggested that geothermal operations in these parts of The Netherlands are unlikely to generate felt seismic events (M>2.0); nevertheless, adequate seismic monitoring techniques are essential to guarantee sustainable and safe use of the Dutch subsurface. Between July and October 2019, Delft University of Technology, Seismotech (Greece), and Gastreatment Services BV installed a passive seismic network to monitor the seismic activity over Nature's Heat geothermal reservoir. The seismic network consists of 30 three-component short-period seismic sensors placed at inter-station distances of approximately 150 m along two crossing lines. A challenge for seismic monitoring systems in urban areas is the high level of background noise. In Kwintsheul, anthropogenic noise dominates the spectrograms at frequencies higher than 2 Hz. Despite these high background-noise levels, a seismic event of ML = 0.0 (duration magnitude Md 0.16) was recorded by all seismometers of the array on July 14, 2019. To understand the relation of the event and improve the safety of the geothermal operation, we are developing a probabilistic monitoring and inversion scheme. This study aims to improve the seismic network's detection and hypocentre-determination capabilities and verifies via template matching if the detected seismic event is repeating over time (possibly at the background noise level). @en

We report on the extraction of deep ocean travel time variations from time-lapse cross-correlations between a hydrophone station and a three-component broadband seismometer. The signals we cross-correlate in this study result from repeated activity by the Monowai seamount, one of the most active submarine volcanoes of the Tonga-Kermadec ridge. In particular, we introduce a specific workflow to exploit repetitive hydroacoustic underwater source activity, which we detail to such an extent that it serves as an example (or “cookbook”). For this reason, we have made the source code publicly available. The workflow proposed in this study (a) overcomes differences in instrument sensitivity and sample rates, (b) involves the selection of eligible cross-correlations based on a source activity criterium as well as slowness analysis, and (c) extracts the travel time variations in distinct frequency bands. In our case, the two frequency bands are 3–6 and 6–12 Hz. We find that the estimated travel time variations in both frequency bands consist of a complex periodic pattern superimposed on a robust linear trend. This linear trend is decreasing, which we attribute to increasing water temperatures along the propagation path of the hydroacoustic signals.

@en

We present an overview of induced seismicity due to subsurface engineering in the Netherlands. Our overview includes events induced by gas extraction, underground gas storage, geothermal heat extraction, salt solution mining and post-mining water ingress. Compared to natural seismicity, induced events are usually small (magnitudes ≤ 4.0). However, due to the soft topsoils in combination with shallow hypocentres, in the Netherlands events exceeding magnitude 1.5–2.0 may be felt by the public. These events can potentially damage houses and infrastructure, and undermine public acceptance. Felt events were induced by gas production in the north of the Netherlands and by post-mining water ingress in the south-east. Notorious examples are the earthquakes induced by gas production from the large Groningen gas field with magnitudes up to 3.6. Here, extensive non-structural damage incurred and public support was revoked. As a consequence, production will be terminated in 2022 leaving approximately 800 billion cubic metres of gas unexploited. The magnitudes of the events observed at underground gas storage, geothermal heat production and salt solution mining projects have so far been very limited (magnitudes ≤ 1.7). However, in the future larger events cannot be excluded. Project- or industry-specific risk governance protocols, extensive gathering of subsurface data and adequate seismic monitoring are therefore essential to allow sustainable use of the Dutch subsurface now and over the decades to come.

@en

Estimating earthquake parameters, including their uncertainty, requires probabilistic sampling or inversion using Bayesian algorithms. One such Bayesian algorithm known to be highly efficient is the Hamiltonian Monte Carlo (HMC) algorithm, and modifying the algorithm with an additional linearization step can further increase this efficiency. However, the modified HMC relies heavily on accurate prior information to effectively sample non-linear earthquake parameters (e.g., hypocenter and origin time). Furthermore, the ability of the modified HMC to estimate non-linear parameters diminishes with respect to the high degree of non-linearity that is inherent to some types of events, such as induced earthquakes. To address this, we adjust the modified HMC to be run in multiple stages, combined with pre-determined initial prior sets. We test this adjustment using synthetic and real data from an induced earthquake event in the Groningen gas field in the Netherlands. We start by obtaining an initial estimate of the prior information and use it to draw multiple initial prior sets. We then run the HMC for each initial prior set in multiple stages where the results from the current stage serve as the prior for the next stage. As the final step, we form the final posterior distributions by selecting results that give the best fit between the observed and modeled data. Within this approach, we estimate ten earthquake parameters those are the six components of a full moment tensor solution, the centroid (three coordinate components), and the earthquake's origin time (including the static time corrections for each recording station). After obtaining the final results, we compare our findings with those of an existing earthquake catalog and several other research results. Given the available fault map of Groningen's subsurface, we found that our results have a higher degree of correlation with respect to the major subsurface faults.@en

We present an efficient probabilistic workflow for the estimation of source parameters of induced seismic events in three-dimensional heterogeneous media. Our workflow exploits a linearized variant of the Hamiltonian Monte Carlo (HMC) algorithm. Compared to traditional Markov chain Monte Carlo (MCMC) algorithms, HMC is highly efficient in sampling high-dimensional model spaces. Through a linearization of the forward problem around the prior mean (i.e., the “best” initial model), this efficiency can be further improved. We show, however, that this linearization leads to a performance in which the output of an HMC chain strongly depends on the quality of the prior, in particular because not all (induced) earthquake model parameters have a linear relationship with the recordings observed at the surface. To mitigate the importance of an accurate prior, we integrate the linearized HMC scheme into a workflow that (i) allows for a weak prior through linearization around various (initial) centroid locations, (ii) is able to converge to the mode containing the model with the (global) minimum misfit by means of an iterative HMC approach, and (iii) uses variance reduction as a criterion to include the output of individual Markov chains in the estimation of the posterior probability. Using a three-dimensional heterogeneous subsurface model of the Groningen gas field, we simulate an induced earthquake to test our workflow. We then demonstrate the virtue of our workflow by estimating the event's centroid (three parameters), moment tensor (six parameters), and the earthquake's origin time. Using the synthetic case, we find that our proposed workflow is able to recover the posterior probability of these source parameters rather well, even when the prior model information is inaccurate, imprecise, or both inaccurate and imprecise.@en

Seismic interferometry (SI) refers to the principle of generating new seismic responses using crosscorrelations of existing wavefield recordings. In this study, we report on the use of a specific interferometric approach, called seismic interferometry by multidimensional deconvolution (SI by MDD), for the purpose of retrieving surface-wave responses. In theory, SI by MDD suffers less from irregularities in the distribution of (passive) sources than conventional SI. Here, we confirm this advantage for the application to surface waves originating from regional earthquakes close to Central Chile. For that purpose, we use the Malargüe seismic array in Argentina. This T-shaped array consists of two perpendicular lines of stations, which makes it rather suitable for the application of SI by MDD. Comparing the responses retrieved through SI by MDD to the responses retrieved using conventional SI, we find that the application of SI by MDD results in surface-wave responses that are both more accurate and more stable than surface-wave responses that are retrieved using conventional SI. That is, our results demonstrate that SI by MDD suffers less from non-uniformly distributed earthquakes and differences in the power spectra of earthquake responses. In addition, we show that SI by MDD mitigates the effect of site amplification on the retrieved surface waves.@en

This paper presents a study on the application of the coda-wave interferometry (CWI) technique using embedded transducers for the monitoring of concrete structural members. The coda-wave interferometry (CWI) has been recognized as a promising technique in detecting small changes in the parameters of materials. To reduce the influence of environment and enable the access to the interested region of the target structure we propose to apply CWI with embedded transducers for the monitoring concrete structures. In this paper, the proposed approach is applied on a 0.5-m cylinder under compression and on a 10-m full-scale beam under bending. The results showed a good correspondence between the CWI-deducted wavespeed change and the estimated strain value in both the small-scale and the large-scale test. Besides, the correlation coefficient of the coda, especially the late coda is sensitive to the formation of cracks even at a distance from the sensor pair. Thus, they can be used as crack indicators of the structure.@en

We present a probabilistic scheme to invert for hypocenters and moment tensors (MTs) of induced seismic events in 2D heterogeneous media. Our scheme is based on a variant of the Hamiltonian Monte Carlo scheme that uses the linearization of a misfit function. This algorithm, however, requires a sufficiently accurate prior, which could be taken from earthquake catalogues. In our case, a detailed subsurface velocity model allows the hypocenter prior to be obtained via a first arrival based algorithm. Fixing the induced event's location to this prior, the MT's prior is subsequently obtained. Having all necessary priors, we then recover the posterior probability density of the induced event's source characteristics, i.e., hypocenter and MT. By comparing the recovered posterior with the directly modelled events, we infer that our scheme is efficient and practical to invert for source characteristics of induced seismic events.

@en

Timing errors are a notorious problem in seismic data acquisition and processing. A technique is presented that allows such timing errors to be recovered in a systematic fashion. The methodology relies on virtual- source responses retrieved through the application of seismic interferometry (SI). In application to recordings of ambient seismic noise, SI involves temporal averaging of time-windowed crosscorrelation measurements. The retrieved interferometric responses are typically dominated by surface waves. Under favorable conditions, these interferometric responses therefore approach the surface-wave part of the medium's Green's function. Additionally, however, its time-reverse is often also retrieved. This implies time-symmetry of the time-averaged receiver-receiver crosscorrelations. In this study, this time-symmetry is exploited: by comparing the arrival time of the direct surface wave at positive time to the arrival time of the direct surface wave at negative time for a large a number of receiver- receiver pairs, relative timing errors can be determined in a weighted least-squared sense. The proposed methodology is validated using synthetic data. The results hold particular promise for large N seismic arrays.@en