Induced seismicity, referring to earthquakes triggered by human activities, has become a pressing concern, particularly in regions with extensive subsurface resource extraction, such as the Groningen gas field in the Netherlands. Decades of natural gas production have altered regional subsurface stress and pressure, reactivating faults and leading to frequent earthquakes, which pose significant structural and societal concerns. To better understand the nucleation and mechanisms behind these earthquakes, source characterization has become a fundamental scientific approach. Accurate characterization not only provides insights into earthquake dynamics but also plays a crucial role in developing subsurface models, optimizing extraction techniques, and helping to mitigate future seismic risks. This research focuses on designing an efficient characterization method by developing a probabilistic workflow for estimating induced earthquake source parameters. Specifically, it integrates advanced probabilistic inversion techniques with geological knowledge to enhance accuracy while reducing computational costs. The proposed workflow leverages a variant of the Hamiltonian Monte Carlo (HMC) algorithms, which offers advantages over a traditional Markov Chain Monte Carlo (MCMC) algorithm, particularly in handling high-dimensional parameter spaces. The workflow estimates key earthquake source parameters, including the centroid, moment tensor, and origin time, even when prior information is limited or uncertain. Initial validation through a synthetic earthquake scenario demonstrates the workflow’s ability to recover source characteristics with high confidence. To further evaluate its validity, the workflow is first tested on a geomechanically simulated synthetic event before being applied to real earthquake data from Groningen. A geomechanical simulation of the 2018 M$_L$ 3.4 Zeerijp event evaluates whether a point source assumption sufficiently represents induced earthquakes compared to finite fault models. The results indicate that, despite simplifications, estimated earthquake parameters remain robust, suggesting that computationally expensive finite fault models may not always be necessary. Building on this, the workflow is then applied to real earthquake data from the 2019 M$_L$ 3.4 event below Westerwijtwerd village, where geological prior knowledge is incorporated to enhance computational efficiency. Using fault-based prior distributions instead of uniform priors significantly reduces computational time while maintaining accuracy, demonstrating the benefits of integrating geological constraints into probabilistic inversion. Expanding its applicability, the workflow is further employed to characterize ten additional earthquakes (M$_L > 2$) in Groningen, where it is combined with rupture directivity analysis to provide insights into rupture dynamics. The results show that normal faulting is the dominant mechanism, aligning with regional geological structures. While rupture directivity effects are minor, they exhibit small correlations with known fault orientations. The study concludes with recommendations for future improvements, such as incorporating S-wave analysis, which could enhance the accuracy of the estimated parameters. Additionally, a more rigorous treatment of uncertainty in both observed and synthetic seismograms is suggested to improve probabilistic estimates. All in all, this research contributes to the advancement of seismic inversion techniques, offering valuable tools for assessing and mitigating the risks associated with induced seismicity. The proposed workflow provides a framework applicable beyond Groningen, aiding in the characterization of induced earthquakes in other regions affected by subsurface activities. ... Read more

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. ... Read more

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. ... Read more

Dynamic geomechanical modeling can generate the seismic wavefield caused by a fault rupture. In dynamic fault-rupture modeling, the source is considered to be finite, with a limited extent both in space and in time. This contrasts with the definition of a point source, which is generally assumed to explain the seismic wavefield caused by an earthquake. Most earlier seismic inversion studies, including those of the induced earthquakes caused by depletion of the Groningen gas field, were performed assuming a point source. Still, finding a point-source reference from the seismic wavefield, even when generated by finite faulting, is important in order to calibrate the geomechanical simulation with field-seismic observations. To this end, we have developed a workflow that links geomechanical forward modeling to seismic moment-tensor inversion. We have tested this workflow for the dynamic rupture considering a realistic 3D layered earth model. At first, we simulate the triggering of dynamic fault slip at the center of a fault plane. Next, we invert the seismograms recorded by receivers located on or near the surface to obtain the full moment-tensor point-source representation and the location of the earthquake. The results of inversion show similar waveforms for both the point source and the finite source. The location of the inverted point source is within 400 m from the center of the slip patch. The double-couple components of the inverted moment tensor also match with the strike and the dip of the fault plane. ... Read more

Geomechanical modelling is generally used to simulate the nucleation of induce d earthquakes in, for instance the Groningen gas field. We apply quasi static simulation to investigate the stress changes from gas production. When a fault reaches a critical state, dynamic simulation provides information on the dynamic rupture during ea rthqu ake nucleation and the resulting wavefield . With the use of geomechanical modelling, it is possible to investigate the effects of the model parameters, e.g., depletion pattern and friction parameters. I n the modelling, the dynamic rupture at a finite fault is simulated both in space and time. The generated seismic wavefield from such a finite source is considered to be more realistic than the resulting wavefield from a point source. T he latter is often assumed in previous studies on the inversion of in duced earthquake data in the Groningen area. To link the wavefield generated by a geomechanically simulated finite source to the field seismic data for an earlier earthquake, we apply the same full moment tensor inversion to the waveform of a finite and of a point source . The inverted moment tensor from the field seismic observation provides a constraint to our geomechanical simulation. This allows us to perform a more realistic simulation of an induced earthquake. ... Read more

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. ... Read more

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. ... Read more

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. ... Read more

The ability to extract information from scattered waves is usually limited to singly scattered energy even if multiple scattering might occur in the medium. As a result, the information in arrival times of higher-order scattered events is underexplored. This information is extracted using fingerprinting theory. This theory has never previously been applied successfully to real measurements, particularly when the medium is dispersive. The theory is used to estimate the arrival times and scattering paths of multiply scattered waves in a thin sheet using an automated scheme in a dispersive medium by applying an additional dispersion compensation method. Estimated times and paths are compared with predictions based on a sequence of straight ray paths for each scattering event given the known scatterer locations. Additionally, numerical modelling is performed to verify the interpretations of the compensated data. Since the source also acts as a scatterer in these experiments, initially, the predictions and the numerical results did not conform to the experimental observations. By reformulating the theory and the processing scheme and adding a source scatterer in the modelling, it is shown that predictions of all observed scattering events are possible with both prediction methods, verifying that the methods are both effective and practically achievable. ... Read more